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a  computer  program  called  PENC02D  has  been  developed  to  analyze  the  non¬ 
normal  Impact  and  penetration  of  any  rigid  axlsywetric  projectile  into  soil 
or  rock.  Horizontal  target  layering  is  allowed  with  accounts  made  for  free- 
surface,  Interface,  and  wake  separation-reattachment  effects.  Expressions  for 
the  compressive  normal  stress  presented  in  U.  S.  Army  Engineer  Waterways  Ex¬ 
periment  Station  Technical  Report  SL-_2fl-a4-fffe  summed  over  the  projectile  sur¬ 
face,  and  the  equations  ofmotion  are  integrated  step-by-step  using  - - 
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20.  ABSTRACT  (Continued) 

a  two-dimensional  Newtonian  integration  scheme. 

Variables  and  format  of  PENC02D's  input  file,  as  well  as  output  options, 
are  given.  To  demonstrate  use  of  PENC02D,  three  sample  problems  are  discussed 
with  setup  and  output  shown. 

A  utility  program  called  MOMENT,  which  calculates  projectile  weight, 
center-of-gravity  location,  and  mass  moments  of  inertia,  is  presented  along 
with  two  example  calculations. 
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CONVERSION  FACTORS,  NON-SI  TO  SI  (METRIC) 
UNITS  OF  MEASUREMENT 


Non-SI  units  of  measurement  used  in  this  report  can  be  converted  to  SI 
(metric)  units  as  follows: 

Multiply _  By  To  Obtain 


g's  (acceleration 
of  gravity) 

9.80665 

metres  per  second  squared 

feet 

0.3048 

metres 

feet  per  second 

0.3048 

metres  per  second 

inches 

0.0254 

metres 

inches  per  second 

0.0254 

metres  per  second 

pounds  (force) 

4.44822 

newtons 

pounds  (force)  per 
square  inch 

6894.757 

pascals 

pounds  (force)- 
seconds  squared 
per  inches  fourth 

10686893.0 

kilograms  per  cubic  metre 

pound  (mass)- 
inches  squared 

0.0002926397 

kilogram-metres  squared 

pounds  (mass)  per 
cubic  foot 

16.01846 

kilograms  per  cubic  metre 

NON-NORMAL  PROJECTILE  PENETRATION  IN  SOIL 
AND  ROCK:  USER'S  GUIDE  FOR  COMPUTER  CODE  PENCQ2D 

CHAPTER  1 
INTRODUCTION 

1.1  BACKGROUND 

In  the  Interest  o£  developing  a  two-dimensional  (2D)  penetration 
code  based  on  easily  obtainable  input,  the  U.  S.  Army  Engineer  Waterways 
Experiment  Station  (WES)  acquired  a  copy  of  Avco  Systems  Division's 
differential-area  force- law  (DAFL)  code  from  Picatinny  Arsenal  in  1975. 
The  code,  in  its  original  form,  uses  a  numerical  predictor-corrector 
scheme  to  solve  the  nonlinear  equations  of  motion  for  a  rigid  projectile 
in  three  dimensions  (3D).  The  differential  force  acting  on  each  surface 
element  is  specified  as  a  function  of  the  local  surface  velocity  and  a 
set  of  empirical  target  parameters.  Avco's  force  law  is  intended  to 
account  for  static  resistance,  surface  friction,  free-surface  chipping, 
target  layer  interfaces,  shock  motion,  and  other  phenomena  in  proportion 
to  the  values  assigned  to  the  various  target  parameters  (Reference  1). 

In  the  hands  of  an  experienced  user,  the  Avco  DAFL  code  is  a 
useful  device  for  extrapolating  laboratory  penetration  test  results  to 
field  predictions.  The  various  target  parameters,  however,  are  only 
vaguely  related  to  standard  engineering  properties  and  are  difficult  or 
impossible  to  measure  directly.  For  the  engineer  at  large,  the  credi¬ 
bility  of  the  code  diminishes  with  the  amount  of  guesswork  needed  for 
the  target  parameters. 

From  1975-1979,  WES  undertook  to  redesign  the  DAFL  code  for  users 
who  are  not  necessarily  penetration  experts.  The  general  structure  of 
the  code  was  changed  in  that  motion  in  3D  was  reduced  to  2D,  making  the 
equations  of  motion  linear  and  replacing  the  predictor-corrector  scheme 
by  simple  Newtonian  integration.  The  main  thrust  of  the  WES  effort, 
however,  was  to  replace  Avco's  differential  force  model  with  a  model 
using  standard  properties  and/or  readily  estimable  empirical  parameters. 


With  the  WES  2D  code,  named  PENC02D,  there  is  less  need  for  guesswork  on 
the  part  of  the  user.  The  major  input  parameters  are  unconfined 
strength  for  hard  targets  (rock,  concrete,  etc.)  and  S -number  for  soil 
(see  Table  1  from  Reference  2) .  The  latter  is  a  penetrability  index 
that  can  usually  be  estimated  within  a  factor  of  2  from  visual  descrip¬ 
tions  of  the  target. 

1.2  SCOPE 

This  report  presents  PENC02D  as  a  black  box;  no  attempt  will  be 
made  to  familiarize  the  reader  with  the  detailed  inner  workings  of  the 
code.  Instead,  input  and  output  variables  (files  and  formats  mainly) 
will  be  defined  and  discussed  (Chapter  2).  To  demonstrate  code  use, 
three  sample  problems  will  be  set  up  and  run  (also  in  Chapter  2) , 
displaying  the  resulting  input  and  output  of  PENC02D.  The  theory  has 
previously  been  presented  (Reference  3);  however,  free  surfaces  above 
and  below  the  penetrating  projectile  are  handled  by  PENC02D  differently 
than  in  Reference  3,  and  are  discussed  in  Appendix  A.  The  2D  Newtonian 
integration  scheme  is  outlined  in  Appendix  B  with  specific  regard  to 
selection  of  time-step  parameters.  A  glossary  of  input  variables  is 
given  in  Appendix  C.  PENC02D  monitors  the  instantaneous  lateral  motion 
throughout  a  given  run  and  will  cut  it  off  completely  if  it  becomes 
insignificant;  guidelines  for  lateral-motion  cutoff  are  discussed  in 
Appendix  D.  Appendix  E  discusses  an  input  parameter  called  ELAS,  used 
to  smooth  the  end  of  the  stress  distribution  and  stabilize  the  numerical 
scheme;  guidelines  are  set  forth  to  help  the  user  select  a  reasonable 
value  for  ELAS.  Appendix  F  introduces  a  utility  program  called  MOMENT 
which  can  be  used  to  calculate  center-of-gravity  (CG)  location,  total 
weight,  and  mass  moments  of  inertia  for  a  penetrator;  these  quantities 
are  required  input  for  PENC02D.  Appendix  G  gives  a  listing  of  PENC02D, 
and  Appendix  H  lists  the  notations  used  herein. 


CHAPTER  2 


USE  OF  PENC02D 


2.1  BACKGROUND 

PENC02D  calculates  the  planar  penetration  of  a  rigid  nonspinning 
axisymmetric  projectile  into  any  target  which  can  be  characterized  using 
either  S-number  (soft  soil-like  materials)  or  density  and  yield  strength 
(hard  rock-like  materials).  Input  includes  projectile  weight  and 
geometry,  target  descriptors,  initial  conditions,  separation  and  free- 
surface  parameters,  integration  time-step  definers,  and  variable 
limiters  to  determine  problem  end  conditions.  Output  is  available  in 
both  printed  and  plotted  format.  Printout  includes  an  echo  print  of  the 
input,  calculated  geometry  table  detailing  each  longitudinal  element's 
shape  and  location,  and  step-by-step  time  history  in  both  the 
projectile-fixed  and  target-fixed  coordinate  systems. 

For  the  most  part,  PENC02D  retains  the  theory  discussed  in  Refer¬ 
ence  3.  The  methodology  for  handling  free  surfaces,  however,  has  been 
modified  somewhat;  a  discussion  of  the  current  treatment  of  layer 
interfaces  and  free  surfaces  appears  in  Appendix  A.  The  original  3D 
version  of  Avco's  DAFL  code  used  a  RUNGE-KUTTA  integration  for  four  time 
steps  to  start  a  predictor-corrector  integration  scheme  used  for  the 
remainder  of  the  penetration  event.  Once  WES's  version  was  cut  to  2D,  a 
simpler  2D  Newtonian  scheme  replaced  the  original  technique  because 
numerical  problems  arose  while  trying  to  integrate  the  current  normal- 
stress  functional  forms  coupled  with  the  often  discontinuous  separation 
and  f  retv-surf  ace  effects.  Details  of  the  current  integration  method  in 
PENC02D,  including  guidelines  for  selecting  related  input  quantities, 
are  given  in  Appendix  B.  It  may  be  useful  for  a  prospective  user  to 
first  examine  Appendices  A  through  E,  which  all  deal  with  PENC02D  input 
quantities  and  give  some  insight  for  selecting  the  more  subjective  input 
parameters . 


2.2  INPUT 


PENC02D  has  one  input  file.  Figure  2.1  is  a  recor--by-record  list 
of  this  file  including  the  required  format  for  each  record.  (Appendix  C 
is  a  glossary  of  PENC02D's  input  quantities  and  contains  the  definitions 
of  all  the  variables  appearing  in  Figure  2.1.)  Each  input  data  record 
is  entirely  alphanumeric,  integer  or  real  (i.e.,  no  mixed-format 
records).  The  first  record  is  the  only  alphanumeric  record  and  contains 
an  identifying  title  for  the  current  problem  in  a  left-justified  3X,A45 
format;  titles  longer  than  45  characters  will  be  truncated  when  printed 
out.  Integer  records  are  read  with  a  1216  format.  Real  records  are 
read  with  either  4E18.6  or  5E14.4  format;  Figure  2.1  gives  details. 

In  addition  to  the  input  guidance  detailed  in  Appendices  A  through 
E,  Figures  2.2  and  2.3  will  assist  the  user  in  the  proper  selection  of 
some  of  the  projectile's  geometric  parameters.  Figure  2.2  details  the 
four  automatic  nose  geometries  available  in  PENC02D.  If  the  projectile 
is  of  one  of  the  types  shown,  the  appropriate  value  of  IOPSHP  is  input 
along  with  values  for  all  of  that  shape's  "required"  parameters  (indi¬ 
cated  in  Figure  2.2)  and  three  other  variables:  DELS(l),  NEM0V2,  and 
ARATIO.  With  this  information,  PENC02D  automatically  divides  the 
projectile  longitudinally  into  several  elements.  A  typical  longitudinal 
element  is  shown  in  Figure  2.3.  DELS(l)  defines  the  thickness  of  the 
first  (nose-tip)  element.  NEM0V2  times  2  defines  the  number  of  circum¬ 
ferential  elements  into  which  each  longitudinal  element  will  be  divided. 
In  Figure  2.3,  DA  represents  one  of  the  surface  differential  areas  that 
results  after  the  circumferential  division  occurs,  with  L  and  W 
representing  the  length  and  width  of  each  DA  on  that  particular  longi¬ 
tudinal  element.  Eliminating  L  and  W  from  the  expression  for  ARATIO 
in  the  figure  and  solving  for  DELS  gives 


DELS(i) 


ARATIO 


COS(THETA(i) 

NEM0V2 


(2.1) 


where  the  subscript  i  refers  to  the  ith  longitudinal  element  from  the 
nose  tip.  This  expression  shows  how  changing  ARATIO  and  NEM0V2  will 
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affect  the  sizes  of  the  longitudinal  elements.  Increasing  ARATIO  and/or 
decreasing  NEM0V2  produces  fewer  longitudinal  elements,  each  having  a 
larger  thickness  and  larger  differential  areas. 

If  the  current  problem  requires  a  nose  shape  different  from  the 
types  shown  in  Figure  2.2,  there  is  provision  for  an  arbitrary  nose 
shape  initiated  by  inputting  IOPSHP  =  0.  The  R-,  THETAD-,  and  ,DELS- 
arrays  must  be  included  at  the  end  of  the  input  file  as  shown  on  sheet  3 
of  Figure  2.1.  NEL  is  the  number  of  longitudinal  elements  to  be  input 
and  is  calculated  automatically  if  IOPSHP  >  0.  The  user  need  not  input 
meaningful  values  of  ARATIO  and  DELS(l)  (sheet  1  of  Figure  2.1)  if  an 
arbitrary  shape  is  to  be  input,  because  neither  quantity  will  be  used. 

Before  concluding  this  section,  it  is  worth  noting  that  PENC02D 
allows  stacking  of  data  sets  in  its  input  file.  In  fact,  execution 
continues  until  an  end-of-file  is  encountered  in  the  input  file.  So,  if 
the  user  has  several  problems  to  run,  either  separate  input  files  can  be 
set  up  and  run  or  one  large  file  with  all  the  data  sets  sequentially 
listed  can  be  run. 

2.3  OUTPUT 

PENC02D  generates  both  printed  and  plotted  output.  The  printout  is 
headed  by  the  title  that  appears  on  the  first  line  of  the  input  file. 

How  much  output  a  user  receives  depends  on  four  groups  of  input  varia¬ 
bles.  The  first  group  specifies  how  many  plots  (NUMPLTS) ,  in  addition 
to  the  two  automatically-generated  plots  of  projectile  shape  and  final 
trajectory,  and  which  plot  types  (ICHOOZ-array)  will  be  output.  The 
second  group  includes  NPRINT,  FREKOT,  and  FREQI  and  defines  the  output 
station  frequency.  The  third  group  is  NSTP,  GAMSTP,  YSTOP,  ZSTOP, 

ALSTP,  TIMEF,  and  VELF  and  encompasses  the  different  ways  to  stop  a 
problem.  The  fourth  group  includes  NUMSTP,  DGAMAD,  and  DALFAD  and 
allows  multiple  executions  with  the  same  data  by  incrementing  initial 
orientation  and/or  attack  angle. 

Plotting  details  are  controlled  with  variables  in  the  first  group. 
If  NUMPLTS  =  0,  the  user  gets  two  automatic  plots.  The  first  is  a  plot 


of  the  resulting  projectile  based  on  the  calculated  R-,  THETAD-,  DELS- 
tables  regardless  of  whether  the  user  inputted  the  tables  or  allowed 
them  to  be  generated  automatically.  The  second  plot  gives  a  trajectory 
of  the  projectile's  C G.  If  NUMPLTS  >  0,  then  an  array  (ICHOOZ)  must  be 
input  having  NUMPLTS  values.  Each  integer  value  of  ICHOOZ  indicates  a 
desired  plot  type  according  to  the  table  shown  under  the  ICHOOZ  listing 
in  Appendix  C.  As  shown  in  the  table,  NUMPLTS  can  be  as  large  as  8 
since  there  are  eight  types  of  additional  plots  available.  If  plot  type 
"8"  is  to  be  selected,  it  must  be  the  last  member  of  the  ICHOOZ-array. 

The  variables  in  group  two  are  important  because  they  determine  how 
many  output  stations  are  to  be  generated.  NPRINT  allows  the  user  to 
select  whether  the  output  locations  will  occur  in  equal  time  intervals 
(NPRINT  =  1)  or  equal  path-length  intervals  (NPRINT  =  2) .  FREKOT  is  the 
size  of  the  interval  in  seconds  if  NPRINT  =  1,  or  in  inches  if  NPRINT  ■ 
2.  FREQI  specifies  how  often,  in  projectile  lengths,  scaled  projectiles 
appear  on  the  automatic  trajectory  plot  discussed  above.  Making  FREQI 
equal  to  a  very  large  number  will  result  in  projectiles  plotted  at  the 
beginning  and  end  of  the  trajectory  only. 

The  variables  in  group  three  give  the  user  several  ways  to  complete 
the  problem.  The  details  of  how  each  variable  works  are  discussed  under 
the  appropriate  listing  in  Appendix  C.  Each  printed  output  station 
shows  the  current  value  of  an  output  variable  called  RETURN.  Table  2.1 
lists  the  possible  values  of  RETURN.  All  output  stations  except  the 
last  one  should  have  RETURN  «  0.  The  last  output  station  always  shows  a 
value  of  RETURN  greater  than  1,  indicating  precisely  why  the  problem 
ended.  The  user  should  check  this  to  make  sure  the  problem  ended  as 
intended . 

The  variables  in  group  four  allow  re-use  of  the  input  data  with  a 
new  projectile  orientation  and/or  attack  angle.  The  first  run  will  set 
y  111  GAMMAD  degrees  and  a  *  ALPHAD  degrees  to  begin  the  problem. 

After  that  problem  is  finished,  y  and  a  are  incremented  so  that  the 
second  run  starts  with  y  *  GAMMAEH-DGAMAD  degrees  and  a  * 

ALPHAD4-DALFAD  degrees.  In  general,  the  ith  run  starts  with 


GAMMAD  +  (i-1)  •  DGAMAD 


(2.2) 


ai  =  ALPHAD  +  (i-1)  •  DALFAD  (2.3) 

where  y^  and  oc^  indicate  initial  values  of  y  and  a  ,  respective¬ 
ly,  for  the  ith  problem  with  the  current  data  set.  NUMSTP  is  the  total 
number  of  runs  to  be  made  (i.e.,  i  takes  on  values  1,2, .. .NUMSTP) .  If 
only  one  run  is  desired  with  a  data  set,  it  is  necessary  to  set  NUMSTP  » 
1.  A  data  set  can  be  skipped  over  by  making  NUMSTP  0. 

The  next  three  sections  will  demonstrate  use  of  PENC02D  by  example 
problems.  Each  problem  will  be  verbally  described,  then  the  input  file 
will  be  set  up,  and  finally  the  resulting  output  will  be  shown. 

2.4  SAMPLE  PROBLEM  1:  SIMULATION  OF 
SANDIA  TEST  R  339-106 

The  first  sample  problem  is  a  simulation  of  a  test  that  was  con¬ 
ducted  by  Sandia  National  Laboratories  in  May  1970  (Reference  4) . 
Although  PENC02D  was  designed  to  handle  2D  problems,  this  problem 
demonstrates  how  it  can  be  used  to  solve  a  one-dimensional  problem 
(i.e.,  normal  impact). 

Figure  2.4  depicts  the  projectile  and  initial  conditions  in  the 

problem.  The  nose  shape  is  an  ogive;  thus,  the  IOPSHP  *  3  automatic 

generating  routine  can  be  used.  From  Figure  2.2  it  is  clear  that  only 

RJ  (the  ogive  radius)  and  R0  (the  aftbody  radius)  need  be  specified, 

in  addition  to  DELS(l),  ARATIO,  and  NEM0V2  to  automatically  generate  the 

longitudinal  elements  (see  Section  2.1).  Because  this  is  a  normal 

impact  problem  with  no  attack  angle  (i.e.,  y  =  180  deg,  a  *  0) ,  it  is 

obvious  that  there  should  be  no  net  lateral  motion.  To  ensure  this,  the 

20  2 

transverse  moment  of  inertia  was  made  very  large  (1.0  x  10  lb-in  )* 
to  avoid  possible  rotation  caused  by  computer  round-off  error.  The 
aftbody  of  the  actual  projectile  used  in  this  test  is  nonflared  and  the 

*  A  table  of  factors  for  converting  non-SI  units  of  measurement  to 
SI  (metric)  units  is  presented  on  page  3. 


motion  is  totally  axial,  so  only  tangential-type  forces  could  occur  on 
the  aftbody.  The  theory  in  PENC02D,  however,  contains  no  tangential 
force  mechanisms  (e.g.,  friction),  and  thus  the  aftbody  can  be  made 
arbitrarily  small.  The  location  of  the  CG  is  also  unimportant  since 
there  is  no  rotation. 

The  target  is  a  three-layer  system  <1  foot  of  3500-psi  concrete 
over  1.5  foot  of  sand  over  a  deep  regime  of  cemented  playa  material). 
Because  there  is  a  hard  layer  over  a  softer  layer,  semifree-surface 
conditions  exist  at  the  concrete/sand  interface  (see  Appendix  A),  and  a 

reasonable  estimate  for  I_  of  the  concrete  should  be  included  in  the 

^  3 

input.  The  density  of  the  concrete  is  estimated  to  be  150  lb/ft 

-3  2  4 

(=0.2246  x  10  lb-s  /in  ).  The  S-numbers  for  the  last  two  layers,  sand 
and  cemented  playa,  are  estimated  to  be  10  and  1.75,  respectively. 

The  projectile  in  the  actual  test  penetrated  about  162  inches. 
NPRINT  was  set  to  2,  so  that  path-length  intervals  are  used  for  output 
incrementing.  In  order  to  obtain  from  100  to  150  output  increments, 
FREKOT  was  set  at  1.25  inches  (i.e.,  162  inches/1.25  inches  per  incre¬ 
ment  =  130  increments).  The  only  plot  of  interest,  besides  the  auto¬ 
matic  shape  and  trajectory  plots,  is  the  graph  of  axial  acceleration 
versus  path  length  (NUMPLTS  =  1,  ICHOOZ(l)  =  8). 

It  is  necessary  to  set  only  one  stop  condition  within  reasonable 
limits,  that  of  final  velocity.  Ideally,  one  would  like  to  run  this 
problem  until  V  =  0.  However,  at  low  velocities  numerical  problems 
arise  that  make  it  impossible  to  run  all  the  way  down  to  zero  velocity. 
The  smallest  tolerable  velocity  is  generally  from  10  to  40  ft/s, 
depending  on  problem  conditions.  For  this  problem,  200  in/s  (=16.67 
ft/s)  will  be  used.  Because  the  last  layer  is  so  hard,  it  is  clear  that 
the  penetrator  would  not  go  much  further  after  a  speed  of  16.67  ft/s 
is  reached  anyway;  thus,  it  is  a  reasonable  end  to  the  problem.  The 
other  stop  conditions  are  set  so  that  they  will  not  interfere  with 
program  execution. 

The  remainder  of  the  input  data  is  either  straightforward  or  was 
achieved  following  guidelines  outlined  in  Appendices  A  through  E.  The 


resulting  input  file,  when  put  together  according  to  the  format  speci¬ 
fied  in  Figure  2.1,  looks  like 


AND1A  SLAB 

8  3 

8 

0 . 250000E 

TEST  R  339-106  (MAY 
195  4  1 

1970) 

1 

60  0  2 

-01  0.200000E 

01 

0 . 100000E-01 

0. 

0. 

0. 

0 . 740000E  02 

0. 

0. 

0 . 400000E 

01 

0. 

0. 

0 . 300000E 

02  0 . 250000E 

02 

0 . 240000E  02 

0. 

0. 

0.834000E 

04 

0 . 125000E  01 

0. 

0 . 200000E 

03  0.100000E 

07 

-0 . 100000E  07 

0. 

0 . 200000E 

00  0 . 100000E- 

02 

0 . 100000E-06 

0 . 100000E  00 

0 . 670000E 

03  0 . 100000E 

21 

0. 

0 . 500000E  01 

0. 

0 . 100000E 

01 

0 . 250000E  01 

0. 

0 . 200000E 

03 

0. 

0 . 180000E  03 

0. 

0.2246E-03  0. 

3500E 

04  0.1200E  02 

0.3500E  03 

1.0000E  01 

0.  0. 

0.3000E  02 

0.1000E  02 

0.1750E  01 

0.  0. 

0.1000E  07 

0.1000E  02 

Figure  2.5  shows  the  first  three  pages,  a  selected  intermediate  page, 
and  the  last  page  of  the  resulting  printed  output.  As  previously 
discussed,  the  first  page  echoes  the  inpyt  while  the  second  page  gives 
the  resulting  longitudinal  sectioning  of  the  projectile,  as  done  by  the 
automatic  generating  routine  for  IOPSHP  -  3.  The  last  element  (I  -  34) 
shows  a  THETAD-value  close  to  -90  deg.  The  automatic  generating  rou¬ 
tines  add  a  thin  (DELS  ■  .001  inch)  element  to  the  tail  end  of  the 
projectile  with  almost  perpendicular  sides  to  close  the  rear  end.  As 
long  as  the  tail  is  closed,  PENC02D  can  handle  situations  in  which  the 
penetrator  tumbles  and/or  flips  over  and  goes  backwards.  Besides,  this 
additional  section  does  not  affect  forward  motion  at  all.  The  column 
labeled  S(I)  is  the  axial  distance  from  the  center  of  each  longitudinal 
element  to  the  CG  of  the  projectile  (positive  if  forward  of  the  CG) . 

The  column  labeled  RCAV(I)  gives  the  local  cylindrical  wake-cavity 
radius  rc  discussed  in  Section  2.4  of  Reference  3  ("Wake  Separation 
and  Reattachment"). 

Sheet  3  of  Figure  2.5  shows  the  first  few  output  stations.  Step  1 
gives  the  initial  conditions  for  the  problem.  Y  and  Z  positions, 
velocities,  forces,  and  g's  refer  to  CG  components  in,Y  and  Z  (target- 
fixed)  directions.  Axial  force  and  g's  refer  to  net  CG  force  and 
acceleration  in  the  z  (projectile-fixed  axial)  direction.  All  the 


quantities  labeled  "PITCH"  refer  to  rotational  motion  about  the  CG: 

ANGLE  -  y,  RATE  -  Y,  FORCE  -  /F„dA,  g's  -  FORCE/m,  MOMENT  -  /ydF  -fz dF 

i  z  y 

(Lquation  B.3),  where  all  integrals  are  taken  over  the  entire  projectile 
surface.  DEL  refers  to  the  time  step  used  just  prior  to  outputting  the 
current  variable  values.  COUNT  indicates  how  many  time  steps  have 
elapsed  since  the  last  output  station. 

Sheet  4  of  Figure  2.5  shows  that  lateral  cutoff  (Appendix  D)  occurs 
between  steps  58  and  59,  because  the  four  rotational  quantities  PITCH 
RATE,  PITCH  FORCE,  PITCH  g's,  and  PITCH  MOMENT  have  all  been  set  to 
zero.  PENC02D  does  not  allow  cutoff  to  occur  until  the  projectile  is 
fully  embedded  in  the  deepest  layer. 

The  last  sheet  of  Figure  2.5  indicates  that,  as  intended,  there 
were  about  130  output  stations.  The  value  of  RETURN  for  the  last  step 
is  7,  indicating  that  the  CG  velocity  has  dipped  below  the  user- 
specified  VELF  (see  Table  2.1),  ending  the  problem  as  desired.  PENC02D 
also  prints  out  final  statistics  consisting  of  total  path  length, 
elapsed  time,  final  velocity,  and  peak  values  of  lateral  and  axial 
accelerations. 

Figure  2.6  gives  the  PENC02D-generated  plots  for  sample  problem  1. 
The  projectile  shape  curve  (sheet  1)  is  provided  for  user  verification. 
The  trajectory  plot  (sheet  2)  indicates  the  path  traveled,  with  projec¬ 
tile  orientation  shown  every  FREQI  projectile  lengths.  In  this  run, 
FREQI  was  input  as  2.5.  The  horizontal  lines  extending  across  the 
trajectory  plot  indicate  layer  interfaces.  Sheet  3  of  Figure  2.6  is  the 
requested  axial  deceleration  versus  path  length  graph  (ICHOOZ  -  8). 

This  particular  curve  is  useful  because  recorded  test  data  are  often 
plotted  this  way. 

2.5  SAMPLE  PROBLEM  2:  BLUNTED 

CONICAL  PROJECTILE  IN  SOIL 

The  second  sample  problem  is  an  oblique  Impact  problem  to  demon¬ 
strate  the  NUMSTP  >  1  option  to  obtain  multiple  runs  with  the  same  data 
set.  The  multiple-run  capability  is  to  be  used  to  search  for  the 


limiting  ricochet  angle  for  a  blunted  conical-nosed  projectile  impacting 
a  soil  target  (S  ■  5)  at  750  ft/s. 

The  penetrator  for  this  problem  is  shown  in  Figure  F.3.  Utility 

program  MOMENT  (Appendix  F)  was  used  to  calculate  weight  (132.55 

pounds),  CG  location  (10.97  inches)  from  nose  tip,  and  transverse  moment 

2 

of  inertia  (3793  lb- in  )  as  listed  in  Figure  F.5.  The  shape  can  be 
program-generated  using  the  IOPSHP  -  4  option  (Figure  2.2)  with  RB  « 

1.0  inch,  RO  ■  3.0  inches,  THETND  *  21.8014  deg,  and  inputting  appro¬ 
priate  values  of  ARATI0(-1.25)  and  NEM0V2(«10). 

Because  this  example  is  concerned  with  projectile  ricochet,  a  wake- 

separation  angle  <J>m^n  *  3  deg  (Reference  3)  is  to  be  used  along  with  a 

f ree-surface  parameter  1^  “  10  (Appendix  A) .  Wake  separation  and  free- 

surface  relief  both  contribute  to  the  rotational  motion  necessary  for 

ricochet;  the  values  chosen  for  4  and  I_  are  typical  for  a  soil 

mm  t 

target. 

An  obliquity  y  “  230  deg  will  be  input  with  the  attack  angle 
a  «  0.  Setting  NUMSTP  ■  3  and  DGAMAD  “  15  deg  produces  three  runs  with 
initial  values  of  y  -  230,  245,  and  260  deg  (Equation  2.2).  Hopefully, 
y  *  230  deg  will  produce  penetration  without  ricocheting,  while  y  - 
245  deg  or  y  -  260  deg  will  produce  ricochet.  If  this  is  the  case,  a 
limiting  ricochet  angle  can  be  interpolated  from  these  results;  or 
another  series  can  be  executed,  starting  with  the  least  oblique  case 
that  did  not  ricochet,  and  setting  DGAMAD  to  a  smaller  increment.  Plot 
types  2,  3,  and  8  are  to  be  selected  (NUMPLTS  -  3,  ICHOOZ  •  2,3,8)  so 
that  the  effect  of  obliquity  on  lateral  and  axial  accelerations  can  be 
examined. 

The  limiting  variable  again  will  be  final  velocity  VELF  (set  - 
200  in/s  for  this  problem),  as  the  purpose  of  the  problem  is  to  run  the 
problem  until  the  projectile  either  stops  or  ricochets.  If  it  ricochets 
(i.e.,  exits  the  soil  halfspace  and  reenters  the  air  layer  above  it), 
PENC02D  will  8 top  the  calculation;  the  printout  will  indicate  this  end 
condition  with  a  value  of  RETURN  -  6  (see  Table  2.1)  in  the  last  output 
station. 


To  select  a  printout  Increment,  it  is  necessary  to  have  at  least  a 
rough  estimate  of  either  the  event's  duration  or  total  path  length.  An 
estimate  of  path  length  can  be  calculated  using  Young's  penetration 
formula  for  normal  impact  (Equation  2  from  Reference  2),  which  is 


Zf  -  .0031  •  S  •  N  •  *  (VQ  -  100)  ,  Vq  >_  200  ft/s  (2.4) 


2 

where  A  is  the  aftbody  cross-sectional  area  («irr0  ).  The  units  for 

2 

this  equation  are  W(pounds),  A(in  ),  VQ(ft/s),  and  Z^(feet).  Sub¬ 
stituting  S  ■  5,  W  «  132.55  pounds,  A  ■  28.27  in^,  and  Vq  *  750  ft/s 
(9000  in/s)  and  estimating  the  nose  performance  coefficient  (see  Table 
II  from  Reference  2)  N  ■  0.6  for  the  blunted  cone.  Equation  2.4  gives 
Zf  -  13.09  feet  (157.1  inches).  If  NPRINT  =  2  and  FREK0T  -  2.0  inches/ 
output  increment,  then  an  estimate  of  the  number  of  output  stations  is 
Z^/FREKOT  *  79  (to  the  nearest  whole  number).  If  the  user  knows  about 
how  many  output  stations  (NQut)  are  desired,  then  Z^/NQut  will  give  an 
estimate  of  FREKOT. 

The  resulting  PENC02D  input  file  looks  like 


BLUNTED  CONE  INTO  S=5  SOIL  (OBLIQUE  IMPACT) 

10  4  198  2  3  3  20  0  2 

2  3  8 

0.  0.125000E  01  0 . 100000E-01  0.150000E  02 

0 . 100000E  01  0.  .  0.  0. 

0.  0 . 300000E  01  0.  0. 

0 . 200000E  02  0.109700E  02  0.500000E  01  0.218014E  02 

0.  0.900000E  04  0.200000E  01  0. 

0 . 360000E  03  0.100000E  07  -0.100000E  07  0. 

0 . 250000E  00  0.100000E-01  0.100000E-07  0.100000E  00 

0 . 132550E  03  0.379300E  04  0.  0.300000E  02 

0.  0 . 100000E  02  0 . 200000E  01 

0 . 300000E  01  0.200000E  03  0.  0.230000E  03 

0.5000E  01  0.  0.  0.1000E  07  0.1000E  02 


Figure  2.7  lists  pages  1,  2,  3,  10,  11,  18,  19,  22,  and  23  from  the 
23-page  printout.  Sheet  2  lists  the  longitudinal  elements  that  result 
from  the  selection  of  ARATIO  "1.25  and  NEMOV2  -  10.  The  best  way  to 
obtain  more  elements  is  to  lower  ARATIO  slightly.  This  is  probably  not 
necessary  because  ARATIO  ■  1.25  means  that  each  DA  is  already  close  to 


being  square  (ARATIO  ■  1.0  would  be  square).  Thus,  although  there  are 
only  23  elements,  this  is  good  enough  for  the  current  problem.  The 
first  element  has  a  value  of  THETA  almost  equal  to  90  deg.  The  reason 
for  this  thin  (DELS  *  .001  inch)  element  is  the  same  as  for  the  rear 
element  being  almost  -90  deg  (explained  for  sample  problem  1,  Section 
2.4).  The  projectile  must  be  a  closed  surface  in  the  direction  of 
travel;  the  bluntness  of  the  nose  requires  a  near  90-deg  element  for 
closure. 

Sheet  3  shows  the  beginning  of  the  first  run  with  step  1  showing 
the  initial  obliquity  y  ■  230  deg.  Note  that  even  at  step  1  there  are 
already  finite  values  for  all  the  forces  and  accelerations.  PENC02D 
orients  the  projectile  so  that  the  nose  tip  sits  initially  at  the  origin 
of  the  target-fixed  coordinate  system,  and  this  results  in  the  initial 
partial  embedment  of  one  corner  of  the  blunted  nose.  To  resolve  this 
problem,  a  user  can  input  a  negative  value  of  ZSH1FT  to  start  the 
calculation  in  the  air  layer  (see  Appendix  C  and  sheet  1  of  Figure  B.2). 

Sheet  4  shows  the  end  of  the  first  run  with  a  value  of  RETURN  ■  7 , 
Indicating  that  the  final  velocity  was  reached.  Thus,  the  projectile 
must  still  be  embedded  in  the  target  as  the  final-position  components 
indicate.  There  were  74  output  stations  for  this  problem  which  is  close 
to  the  prerun  estimate  of  79. 

Sheet  5  initiates  the  second  run  with  the  same  data  set  with 
y  -  243  deg.  The  initial  axial  and  pitch  accelerations  are  both  larger 
in  magnitude  than  the  case  for  y  ■  230  deg  because  more  of  the  nose  is 
initially  embedded  as  the  angle  of  entry  approaches  the  horizontal. 

Sheet  6  ends  the  second  run,  again  indicating  that  the  final  velocity 
was  reached  (RETURN  *  7)  with  the  projectile  still  buried  in  the  target. 

Sheet  7  starts  the  last  run  of  the  sequence  with  y  ■  260  deg.  The 
projectile  is  nearly  horizontal  (y  *  270  is  exactly  horizontal).  As 
expected,  the  initial  axial  and  pitch  accelerations  are  higher  than  the 
first  two  runs  because  a  larger  section  of  the  nose  is  initially 
embedded. 

The  most  interesting  occurrence  appears  in  step  37  (sheet  8) ,  when 
all  the  forces  and  accelerations  suddenly  go  to  zero.  The  position 


coordinates  indicate  that  the  projectile  has  ricocheted  and  reentered 
the  air,  but  the  calculation  continues  until  the  automatic  check  finally 
shuts  down  activities  in  step  46  (sheet  9) .  This  is  verified  with  the 
value  of  RETURN  -  6. 

Figure  2.8  shows  the  requested  plots  for  this  entire  problem. 

Sheet  1  shows  the  shape  graph,  verifying  the  outer  surface  geometry  of 
the  projectile  portrayed  in  Figure  F.3.  Sheets  2  through  5  represent 
output  from  the  y  =  230  deg  run,  sheets  6  through  9  the  y  *  245  deg 
run,  and  sheets  10  through  13  the  y  -  260  deg  run. 

Comparing  the  three  trajectory  plots  (sheets  2,  6,  and  10),  the 
increased  initial  embedment  of  the  projectile  with  increasing  initial  y 
is  clearly  evident,  along  with  the  ricochet  for  y  ■  260  deg.  Based  on 
these  three  runs,  the  limiting  ricochet  angle  (the  smallest  y  for 
ricochet)  must  be  between  245  and  26^  deg.  The  foregoing  procedure 
could  be  repeated,  starting  with  y  ■  245  deg  and  Incrementing  by  5  deg, 
to  more  closely  define  the  limiting  angle. 

Sheets  3,  7,  and  11  are  the  three  axial  deceleration-time  histories 
(1CH00Z  *  2)  requested.  Comparison  of  these  three  graphs  shows  a 
decreasing  peak  axial  deceleration  as  y  Increases.  Similarly,  sheets 
4,  8,  and  12  give  the  three  lateral  acceleration-time  histories 
(1CH00Z  ■  3)  showing  the  negative  extreme  growing  with  increasing 
initial  obliquity.  As  y  gets  closer  to  horizontal,  the  projectile 
encounters  greater  lateral  forces  from  below  due  to  the  top-side  relief 
from  the  free  surface.  Because  ^  >  0,  there  is  a  cavity  around  the 

bottom  side  of  the  aftbody,  and  the  unbalanced  loads  can  easily  cause 
rotation  to  occur.  With  the  increased  angle  of  attack,  the  axial  forces 
must  fall  off  more  rapidly  because  these  lateral  loads  are  causing 
significant  kinetic  energy  loss  from  the  projectile. 

2.6  SAMPLE  PROBLEM  3:  SIMULATION  OF 

SANDSTONE  3-DEG  REVERSE  BALLISTIC  TEST 

In  1977,  Eandia  National  Laboratories  conducted  four  reverse 
ballistic  tests  (RBT)  into  sandstone  for  the  Defense  Nuclear  Agency 
(DNA).  Reference  5  contains  a  complete  description  of  the  test  series. 


The  last  sample  problem  is  a  simulation  of  two  of  the  tests  (Tests  2  and 
3)  which  were  done  at  a  3-deg  angle  of  attack.  Figure  2.9  indicates  the 
nominal  impact  conditions  for  the  tests:  Vq  *■  1500  ft/s  (18,000  in/s), 

Y  *  183  deg,  and  a  =  3  deg.  Figure  2.10  shows  the  surface  geometry  of 
the  actual  test  projectile.  This  sample  calculation  omits  the  aftbody 
flare  (from  diameter  *  1.702  to  1.900  inches).  The  automatic  generating 
routines  can  handle  a  flared  aftbody  if  the  diameters  behind  and  in 
front  of  the  flare  section  are  the  same.  The  RBT  projectile,  however, 
flares  out  to  the  new  diameter  (1.900  inches)  and  then  holds  that 
diameter  constant.  In  the  latter  case,  the  user  must  input  IOPSHP  -  0 
and  input  the  entire  R-,  THETAD-,  and  DELS-arrays.  For  this  example, 
IOPSHP  *  3  will  be  used  to  show  the  automatic  generating  routine's 
ability  to  place  the  conical  tip  on  the  ogival  nose  (see  Figure  2.2  for 
IOPSHP  -  3).  The  flared  aftbody  will  be  replaced  by  a  straight  aftbody 
having  a  diameter  of  1.702  inches. 

Because  of  the  dimensions  of  the  test  target,  only  the  first  1.2  ms 
of  the  experiments  were  considered  valid.  Thus,  the  limiting  variable 
for  this  problem  will  be  TIMEF.  If  NPRINT  ■  1  (to  set  output  increment 
units  to  time)  and  FREKOT  *  0.00001  second,  then  there  should  be  122 
(«1+(TIMEF-TIMEI) /FREKOT)  output  stations  if  TIMEI  -  0  and  TIMEF  - 
0.00121  second. 

Of  interest  in  this  problem  are  axial  deceleration-,  lateral 
acceleration-,  and  total  CG  moment-time  histories,  because  the  projec¬ 
tiles  in  this  test  series  had  instrumentation  to  measure  axial  and 
lateral  accelerations  as  well  as  axial  strains  at  several  locations.  To 
obtain  these  plots,  NUMPLTS  has  been  set  to  3  with  the  ICHOOZ-array 
equal  to  2,  3,  and  4. 

Because  it  is  necessary  to  obtain  more  detail  in  this  problem  than 
in  most  penetration  problems  (i.e.,  FREKOT  is  small),  the  time-step 
parameters  FRAD,  FANG,  NUMN0S  (Appendix  B)  are  much  smaller  than  in  the 
first  two  sample  problems.  They  must  be  set  small  enough  so  that  one  or 
more  Integration  steps  occur  between  each  printout  step.  In  this 
problem  it  was  decided  to  set  the  time-step  parameters  so  that  7  to  11 
integration  steps  would  occur  during  each  printout  increment.  For  most 


problems,  the  time-step  parameters  can  be  set  much  larger  and  still  give 
a  converged  solution. 

In  any  event,  with  all  the  above  information  taken  into  account, 
the  following  input  file  was  used: 


SANDSTONE  RBT  SLED  TEST . . . ALPHA=3  DEG 

8  3  195  2  3  1  150  0  1 

2  3  4 

0 . 250 000E-01  0 . 200000E  01  0.200000E-01  0. 

0 . 176000E  00  0.  0 . 102120E  02  0. 

0.  0.851000E  00  0.  0. 

0 . 181500E  02  0.871000E  01  0.379700E  01  0. 

0.  0 . 180000E  05  0 . 100000E-04  0. 

0 . 100000E  07  0 . 100000E  07  -0.100000E  07  0. 

0 . 200000E-01  0.300000E-03  0.200000E-07  0.200000E-04 

0 . 948000E  01  0 . 221900E  03  0.300000E  01  0.500000E  02 

0.  0 . 121000E-02  0 .200000E  01 

0.  0 . 240000E  03  0.  0.183000E  03 

0.  0 . 1947E-03  0.3400E  04  0.1000E  07  0.3500E  03 


Figure  2.11  lists  the  first  three  pages  and  the  last  page  of  the 
resulting  printout.  Sheet  1  verifies  the  input  discussed  above.  Sheet 
2  reflects  the  longitudinal  sectioning  resulting  from  the  IOPSHP  -  3 
shape  routine.  The  conical  nose  tip  is  evident,  as  the  first  six 
elements  all  show  a  constant  surface  angle  THETAD  of  50.3  deg.  Sheet  3 
is  the  beginning  of  the  problem.  Step  1  shows  a  small  initial  Y- 

_3 

velocity  (.5  x  10  in/s)  when  it  should  be  identically  equal  to  zero, 
indicating  the  presence  of  computer  round-off  error.  The  last  sheet  is 
the  end  of  the  problem  (time  equal  to  TIMEF).  RETURN  ■  3  in  step  122 
verifies  that  TIMEF  was  the  source  of  problem  termination. 

Figure  2.12  exhibits  the  requested  plots  for  this  problem.  Sheets 
1  and  2  automatically  verify  the  projectile  shape  and  path  traveled. 
Sheets  3,  4,  and  5  show  the  axial  deceleration-,  lateral  acceleration-, 
and  CG  moment-time  histories,  respectively. 


Table  2.1  Values  of  RETURN  in  PENC02D  printout. 


EXPLANATION 


Code  executing  normally;  no  stop  conditions  encountered 
yet. 

Maximum  number  of  printout  increments  (NSTP)  has  been 
reached;  calculation  stops. 

Current  time  t  has  exceeded  TIMEF;  calculation  stops. 

Absolute  value  of  Y  (horizontal  coordinate  of  CG)  has 
exceeded  YSTOP;  calculation  stops. 

Value  of  Z  (vertical  coordinate  of  CG)  has  become  less 
than  ZSTOP;  calculation  stops. 

Projectile  has  exited  the  earth  target  and  passed  up 
into  the  air  (layer  1);  calculation  stops. 

CG  velocity  V  has  becbme  less  than  VELF;  calculation 
stops. 

Orientation  angle  y  has  exceeded  GAMSTP;  calculation 
stops. 


Absolute  value  of  angle  of  attack  o  has  exceeded 
ALSTP;  calculation  stops. 


TITLE 


3X.A4S 


IOPSHP 

NSTP 

I4LAY 

NUMPLTS 

NUMSTP 

NUMNOS 

IREDUC 

HPRIHT 

16 

16 

16 

16 

16 

16 

16 

16 

16 

'  OELS(l) 

ARATIO 

ELAS 

E18.6 

El  8. 6 

E18.6 

DGAMAO 


E18.6 


EJ 

RJ 

E18.6 

E18.6 

RO 

XFA 

XFF 

E18.6 

E18.6 

E18.6 

XCG 

SN 

THETND 

E18.6 

E18.6 

E18.6 

VEL 

FREKOT 

DALFAD 

E18.6 

E18.6 

E18.6 

THETFD 


E18.6 


Figure  2.1  Data  format  for  input  file  S  (Sheet  1  of  3). 
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YSTOP 

ZSTOP 

ZSHIFT 

E18.6 

E18.6 

E18.6 

FANG 

DTMIN 

DTMAX 

E18.6 

E18.6 

E18.6 

XICG 

ALPHAO 

ALSTP 

E18.6 

El  8. 6 

E18.6 

TIMEF 

FREQI 

E18.6 

E18.6 

VELF 

Mil 

E18.6 

E18.6 

GAMMAD 


E18.6 


Note:  Include  this  record  for  each  target  layer  I. 

2  <  I  <  NLAY  (Layer  1  Is  air  and  need  not  be  specified) 


E14.4 


End  here  If  IOPSHP  f  0. 


DENSITY(I) 

YIELD(I) 

E14.4 

E14.4 

Note:  Include  enough  records  to  exhaust  entire  THETAD-array. 


^  THETAD(l) 

THETAD(2) 

•  •  • 

THETAO(NEL) 

E18.6 

£18.6 

E18.6 

E18.6 

Note:  Include  enough  records  to  exhaust  entire  R-array. 

^  R(D 

R(2) 

•  •  • 

R(NEL) 

E18.6 

E18.6 

E18.6 

El  8. 6 

Note:  Include  enough  records  to  exhaust  entire  OElS-array. 


Value  of 
IOPSHP 


Shapes  Showing 
Required  Input 


Remarks 


1 


SN  (nose  length)  may  be  Input  Instead 
of  THETND;  THETND  will  then  be  calcu¬ 
lated;  If  THETND  is  Input,  SN  will  be 
calculated. 


2 


AH  Input  shown  is  required. 


3 


(a)  If  a  conical  tip  Is  desired,  all 
Input  shown  Is  required. 

(b)  If  no  conical  tip  Is  desired,  set 
RB  =  0;  input  either  SN  or  RJ  and 
the  other  will  be  calculated;  all 
other  Input  shown  Is  required. 


4 


All  input  shown  is  required;  SN  is 
calculated;  DELS(1)  Is  not  required 
input. 


Figure  2.2  PENC02D  projectile  shape  options. 


XLP  =  77.5  IN 

RO  =  4.0  IN 

WEIGHT  =  670  LB 

RJ  =  74.0  IN 

(OGIVE  SHAPE) 


12- IN  CONCRETE 
(3500  PSI) 


18- IN  SAND 


CEMENTED  DRY  LAKE  PLAYA 
(300  BLOWS/FT) 


Figure  2.4  Initial  conditions  for  sample  problem  1. 
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Figure  2.7  (Sheet  3  of  9) . 
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Figure  2.7  (Sheet  7  of  9). 


! 


1 


i 


OHO 

oo 

OHO 

OO 

OHO 

OO 

Oho 

oo 

♦MO 

OO 

OO 

MUA 

wmh 

IUUAA 

MiMA 

MMA 

MJMMA 

jSSS 

o«S32 

rsMCMA 

INO<NOCO 

OOtfl 

INOOHN 

OsTI/l 

AlOONN 

OHIA 

MOOON 

«r 

sr 

sr 

ST5 

~S£! 

eo 

M  MOO 

M  MOO 

M  MOO 

M  »  OO 

M  MOO 

So® 


ZZZZtC  HtCXI  aCZZZZ 

uuuuo  uooou  uuuuu 


F*JW^  P+fifW*4  iftlAl«KO  H^mo  rfttAf^O 


♦<rt 

•»<Tr4>flN 

*-n/ms»£k 

23ESS 


ftlOy*jrv  WVOHCMH  »HNNIA  ^r»H|v«  <CM> 

HOHtVO  «<MHK)  IAO‘<MH  DOHfSf  HOAKKO  NW 

OTMFSOK  *>HKM/VH  ww  KI^O«rt  ^ 

NNAtH<0  (AOIO^IA  «MKW4H  fsAI 

NOCJTVO  «OHH<D  N U)N  OH 

KIM^<0  KXS^HWH  ««iAHN  ^MHHH  IMlNrM  ON 


<K4> 

OKSi 

MO 

ON 

♦M 


S3 


II  II  N  M  H 

5k  cn 

Comi  o 


M  II I . II  II  II  II  II  II  II  II  II  II  II  II  II  II  II  II  II  II  II  n  N  N  II  N  II  II  II II  M  II II II  II  II  N  II  N  N  N 


«UK«U 

8S£i2 

M  •  IX 

NNNN< 


s~  an  §£ 

KJU  <D  HJHI 
MOW  HOU 

VUtttfU 


5£  »  5f  »  g: 


OuiO-  < 
kl>in<dm 


I  I  I  IX  I  I  I  IX  I  I  MX 

■julJlM  knjtulku# 


vuotfu  vtJaMJ  wjkmj  tfuaiAJ 
OUJO*  <  OUJO-  <  Otuo  <  QUJO-  < 
£>U-0*H  IL>lft.dH  &XL0H  OXkttH 


I  I  I  IX  I  I  I  IX 


ooooo 

II  11  II  II  M 

M  M  II  II  II 

ot  in 

S£  an 

h&  • 

ais  • 

£^&3 

w 

o 


M  I  IX  m  ix  i  i  i  ix  i  i  i  ix 

***■*■>■  -  m  m  iTi  i  m  >■»■»■>•  — 

npw^p^ 


NOtflW I  NSlAKVt  N^ONt  N«ON«  NJ»W  «ft  AW 


<S<6NOO  NIW«^)  N^NOIft  NOON®  <OOOHN  AlMXM  UW 

OAnNSH  KMOOlAN  ONHON  NrtON^  <0O®«O  O® 

HIANO®  ®0®0®  WAKOKHO  OHOOift  ONNO®  OOONO  ®W 

N«®lrtH  <4r»<MtfK«  OtfKVJtfWM  (NsrKKlfl  WK*®K\*M  4)K>®ON  KN1 


oo«ho  ari^o  awwhh 

ONNNH  IfINAM®  IflNHHIft 


-SOON®  ®OOHN  «0®NN  eo 
WMWN  IflNflOH  IANNAN  ON 


32 

MJMI 

3* 

AHA 


S3 

an 

WA 

oo 

♦k 


32 

Mil 

NO 

#t 


I  I  IX  I  l  I  I 


i*i 

§ti  g 

^  Ml 

£  I 

Ml 

£:  S 

MM  O 

Ml 

St  a 

MM  O 

„  Ml 

St  g 

MM  O 

Ui 

st  g 

MM  O 

Ml 

St  g 

MM  O 

KOMI  IN 
MOO  _  . 

353  “■ 

KOMI  IN 

t8o  “■ 

KOMI  IN 
MOO  _ 

KOMI  IN 
MOU 

KOMI  IN 
MOO 

S£b2 

£3 

a.>!NQM 

lAWfttfU 

22£»2 

kAHttirU 

S&H 

bsk»2 

ix  i  i i ix  i  i  i  ix 


IX  l  l  l  IX  I  l  I  »X  i  I  I  ix  i  I  i  IX  1111 


«  c 


a 


i 

K 

M 


M 
IN 

e 

M  M 


W  Ml 

S  I 

3  *  p;  • 


2 

_  IN 

O  O 


w  H  it  Si 

S  P  S  P  s» 


*  ^  %*.  2*  P;  ?  2®  *  *  9* 

N 

!c  SH  S:  S'  Si  *  Si  2  Si  £  Si  2  Si 

W  K  N  H  N  K  M  N  H  h  N  h  H  h 


42 


Figure  2.7  (Sheet  8  of  9). 
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APPENDIX  A 


INTERFACE  AND  FREE-SURFACE  EFFECTS 

The  penetration  resistance  near  a  layer  interface  or  free  surface 
will  usually  differ  from  that  of  a  purely  homogeneous  medium.  The 
PENC02D  code  acknowledges  the  presence  of  such  boundaries  in  a  manner 
that  is  crude  but  qualitatively  correct. 

Each  target  layer  is  assigned  a  parameter  IT  which  loosely 
corresponds  to  the  rigidity  index  (Reference  6)  in  tension.  For  soft 
materials  (soil)  the  recommended  value  is  IT  =  10,  and  for  hard  mate¬ 
rials  (rock  and  concrete)  the  recommended  value  is  1^  =  350.  In  any 
case,  there  corresponds  to  IT  a  damage  radius  rQ  defined  at  each 
point  on  the  projectile  surface 

rD  ■  rp<Vl/3  <*•» 

The  quantity  rp  represents  the  local  cylindrical  radius  of  the  projec¬ 
tile.  At  each  point  on  the  projectile,  a  damage  vector  is  defined  such 
that 

-+  -y  ,  1/3 

rD  - n  yy  <a*2> 

where  n  is  a  unit  vector  normal  to  the  projectile  surface.  The  target 
properties  to  be  used  in  defining  the  normal  stress  acting  on  any 
differential  area  on  the  surface  of  the  projectile  depend  upon  the 
particular  layer  in  which  the  tip  of  the  damage  vector  lies.  In  Figure 
A. 1  the  points  A',  B',  and  C'  are  the  tips  of  damage  vectors  emanating 
from  differential  areas  centered  at  points  A,  B,  and  C,  respectively. 

If  the  tip  is  in  the  same  layer  occupied  by  the  center  of  the  projectile 
surface-area  element  (e.g.,  point  B'  in  Figure  A.l),  the  stress  is 
calculated  in  the  usual  manner  for  that  layer.  But  if  the  tip  falls  in 
an  adjacent  layer  (e.g.,  points  A'  and  C*  in  Figure  A.l),  the  stress 


I 


calculation  is  made  using  the  properties  of  the  softer  of  the  two 
layers . 

The  foregoing  scenario  represents  an  attempt  to  incorporate  the 
effects  of  chipping  and  spallation  into  PENC02D.  For  example,  when  a 
projectile  penetrates  a  soft  layer  followed  by  a  hard  layer,  the  latter 
has  no  effect  until  contact  is  actually  made.  On  the  other  hand,  in  the 
case  of  hard  followed  by  soft,  the  soft  layer  begins  to  have  an  effect 
as  the  interface  is  approached.  In  general,  whenever  a  projectile 
surface  element  is  near  an  interface  or  free  surface,  there  may  be  a 
reduction  in  the  resisting  stress  if  the  adjacent  layer  is  softer  than 
the  immediate  layer. 


» 


I 


LAYER 


APPENDIX  B 


PENC02D'S  NEWTONIAN  INTEGRATION 
WITH  GUIDELINES  TO  RELATED  INPUT 


PENC02D  uses  a  two-dimensional  (2D)  Newtonian  scheme  to  integrate 
the  following  equations  of  motion: 


«  mY  * 

F  cos 

y  -  F  sin  y 

(B.l) 

y 

z 

=  jdZ  ■ 

F  sin 

y  +  F  cos  y 

(B.2) 

y 

z 

-  IY  - 

lyuz  ■ 

-  /zdF 

y 

(B.3) 

Equations  B.l  through  B.3  are  Equations  2.9  through  2.11  from  Reference 
3,  rewritten  using  PENC02D's  coordinate  system.  Figure  B.l  shows  the 
difference  between  the  two  coordinate  systems.  The  integrals  in  Equa¬ 
tion  B.3  are  evaluated  over  the  projectile's  entire  surface  area. 

Figure  B.2  is  a  flowchart  of  the  solution  procedure  used  by  PENC02D. 

This  integration  begins  by  initializing  Y  ,  Z  ,  and  y  (the  two 
spatial  coordinates  of  the  projectile  center  of  gravity  (CG)  and  the 
penetrator  orientation)  and  their  respective  velocities,  as  well  as  time 
t  .  The  current  speed  is  calculated  and  used  to  compute  a  time  step  At 
according  to  the  flowchart  in  Figure  B.3.  Subroutine  FORCES  is  called 
to  sum  the  stress  components  over  the  projectile's  surface  to  obtain 
total  force  components  Fy  and  Fz  and  total  CG  moment  M  .  The 
changes  in  the  translational  and  rotational  velocities  can  then  be 
calculated  for  the  current  time  interval.  New  velocity  components,  as 
well  as  average  velocities  over  the  interval,  are  then  determined  so 
that  the  displacement  component  changes  over  the  current  interval  can  be 
computed.  Once  these  components,  along  with  the  time,  are  updated  and 
all  stop  conditions  checked,  the  integration  loop  cycles  to  the  next 
time  interval.  The  problem  continues  until  one  of  the  inputted  stop 
conditions  is  met. 


B1 


The  following  parameters  are  responsible  for  adjusting  the  current 
time  step  to  be  used  in  the  problem  of  interest:  NUMNOS,  FRAD,  FANG, 
DTMIN,  DTMAX,  and  1REDUC.  Initially,  NUMNOS  is  used  to  embed  the  pene- 
trator  nose.  The  nose  is  divided  into  NUMNOS  equal-thickness  sections, 
and  each  of  the  first  NUMNOS  time  steps  embeds  one  of  these  sections. 

If  S(t)  is  the  path  length  traveled  by  the  projectile  up  to  time  t  , 
and  is  the  projectile  nose  length,  then 


S(ti)  -  s<ti_1)  =*  jjuifljos 


(B.4) 


during  the  nose-embedding  process.  For  a  small  interval  At  , 


ViAti  *  S(ti)  -  S(t±^) 


(B.  5) 


where  is  the  velocity  at  the  beginning  of  the  current  interval. 

The  first  part  of  the  penetration  event  is  thus  characterized  by  a  time 
step  given  by 


At 


i  V±  •  NUMNOS 


(nose  embedment)  (B.6) 


Generally,  10  to  50  is  a  reasonable  range  of  values  for  NUMNOS,  depend¬ 
ing  on  target  material,  impact  velocity,  and  nose  length. 

The  time  step  used  in  the  next  portion  of  the  penetration  event  is 
characterized  using  FRAD  and  FANG.  FRAD  is  used  to  calculate  a  time 
step  At  that  will  translate  the  penetrator  by  a  distance  equal  to 
FRAD  projectile  radii,  while  FANG  is  used  to  calculate  a  At  that  will 
rotate  the  penetrator  by  an  angle  equal  to  FANG  radians.  PENC02D  uses 
the  smaller  of  the  two  time  steps,  allowing  both  rotational  and  transla¬ 
tional  limitations  to  be  set.  Some  problems  require  that  the  rotation 
in  a  given  time  interval  be  limited  (e.g.,  problems  involving  high 
angular  accelerations)  while  other  problems  rely  on  translational  motion 


co  determine  time-step  size.  The  time  interval  can  be  summarized  as 
follows  after  nose  embedment 


(B.7) 


where  rQ  is  projectile  radius  and  and  are  the  current  CG  and 

angular  velocities,  respectively.  How  big  should  FRAD  and  FANG  be?  In 

multilayer  problems,  take  the  thinnest  layer  to  be  penetrated  («T  .  ), 

min 

estimate  the  number  of  time  steps  desired  in  this  layer  (-N^),  and 
calculate  FRAD  as 


FRAD 


r 

o 


(multilayer  problem)  (B.8) 


In  homogeneous  problems,  decide  how  many  steps  are  required  to  accurate¬ 
ly  integrate  the  equation  of  motion;  generally,  NQ  -  100  to  500  are 
required  depending  on  target  material,  projectile  size,  and  impact 
conditions.  Let  Tmin  be  an  estimate  of  the  total  path  length  and  use 
Equation  B.8.  This  will  give  a  starting  value  for  FRAD.  In  all  prob¬ 
lems,  two  or  more  values  of  FRAD  should  be  tried  to  check  for  conver¬ 
gence.  If  a  problem  is  not  dominated  by  rotation,  FANG  can  be  made 
large  enough  so  that  it  will  not  affect  Che  time  interval  calculation. 
Most  of  the  current  experience  is  with  problems  of  relatively  low 
angular  motion.  To  select  a  value  for  FANG,  estimate  the  maximum 
orientation  change  |AyjJ  that  can  be  tolerated  during  one  time  inter¬ 
val  for  the  problem  at  hand.  FANG  is  this  value  in  radians. 

When  the  velocity  gets  close  to  the  final  velocity,  PENC02D  fixes 
the  time  step  at  a  constant  value  to  finish  the  calculation.  Specifi¬ 
cally,  it  uses  the  At  of  the  increment  just  before  gets  within 

2AV^  of  the  final  velocity  (VELF).  IREDUC  is  used  when  the  standard 
termination  procedure  is  inappropriate,  such  as  when  significant  lateral 
motion  still  remains  at  the  end  of  a  problem.  If  IREDUC  >  0  ,  then 


IREDUC  reductions,  each  by  a  factor  of  3,  will  be  made  in  both  FRAD  and 
FANG  when  the  current  CG  velocity  has  fallen  to  certain  fractions 

of  the  impact  velocity  Vq  ,  as  given  by  Table  B.l.  If  lateral  motion 
has  been  cut  off  (Appendix  D)  before  <_0.2Wq  ,  then  these  reductions 
do  not  occur. 

Throughout  the  entire  penetration  event,  the  time  step  will  never 
be  less  than  DTMIN  nor  greater  than  DTMAX,  regardless  of  the  values  of 
FRAD,  FANG,  etc.  If  for  some  reason  a  constant  time  step  At  is 
desired  for  a  particular  problem,  just  set 

DTMIN  ®  DTMAX  =  At  (constant  time  step)  (B.9) 

Otherwise,  DTMIN  and  DTMAX  should  be  used  to  keep  the  time  interval  from 
getting  too  big  or  too  small. 


Table  B.l  Velocities  at  which  integration  parameters  FRAD 
and  FANG  are  revised  using  time-step  reduction. 


Reduction 

Level 

Velocity 

New  FRAD 

New  FANG 

First 

0.2*V 

FRAD  /3* 

FANG  / 3** 

o 

o 

0 

Second 

0.1*V 

FRAD  /9 

FANG  /9 

o 

o 

o 

Third 

0. 05*V 

FRAD  /27 

FANG  /27 

o 

o 

o 

Fourth 

0.025*V 

FRAD  /81 

FANG  /81 

o 

o 

o 

etc. 

etc. 

etc. 

etc. 
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PENC02D 

Variable  Explanation 


PENC02D 


Reference  3 


Projectile-fixed  lateral  direction, 
positive  on  clockwise-side  of  z-axis 


Projectile-fixed  axial  direction, 
positive  forward  of  CG 


Target-fixed  horizontal  direction 
positive  to  the  right 


Target-fixed  vertical  direction, 
positive  upward 


Projectile  orientation  angle, 
measured  from  positive  z-axis 
clockwise  to  upward  vertical 


Y 


Z 


X 


-Z 


Y 


0  +  ir 


PENC02D 


Reference  3 


INITIALIZATION 


t  -  to 

Y  =  Y0  =  XCG  •  sin(Y0) 

1  *  zo  ■  XCG  *  cos(Y0)  -  Zs 

Y  =  Y0  =  -V0  •  sin(Y0  -  a0) 
^  “  ^0  “  Yq  cos(y0  -  c«q) 

Y  =  Yo 


Figure  B.2  Flowchart  for  PENC02D's  Newtonian 
integration  scheme  (Sheet  1  of  3). 
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APPENDIX  C 


ALPHAD 

ALSTP 

ARATIO 

DALFAD 

DELS (I) 
DENSITY (J) 

DGAMAD 

DTMAX 

DTMIN 

ELAS 

EJ 

FANG 

FRAD 

FREKOT 

FREQI 

GAMMAD 


GLOSSARY  OF  INPUT  VARIABLES 


Initial  angle  of  attack,  deg 

Problem  ends  if  absolute  value  of  angle  of  attack  a 
ever  exceeds  ALSTP,  deg 

A  differential  area's  length-to-width  ratio;  used  by 
the  geometry-generating  routines  if  IOPSHP  >  0, 
dimensionless 

Used  if  NUMSTP  >  1,  to  increment  ALPHAD  for  subse¬ 
quent  run(s)  using  current  data,  deg 

Thickness  of  the  ith  longitudinal  element,  in 

Target  mass  density  of  jth  layer;  used  for  rock-like 

lb-s2 
targets,  - r~ 

in 

Used  if  NUMSTP  >  1,  to  increment  GAMMAD  for  subse¬ 
quent  run(s)  using  current  data,  deg 

Maximum  allowable  time  step,  sec 

Minimum  allowable  time  step,  sec 

Normal-stress  smoothing  parameter  (Appendix  E) ,  rad 
Secant  ogive  offset,  in 

Rotational  time-step  parameter  (Appendix  B) ,  rad 

Translational  time-step  parameter  (Appendix  B) , 
radii 

Print  output  increment,  seconds  if  NPRINT  *  1, 
inches  if  NPRINT  *  2 

Distance  along  trajectory  path  between  successive 
scaled  projectiles  on  automatically-generated 
trajectory  plot,  projectile  lengths 

Initial  projectile  orientation  measured  from  the 
line  connecting  the  center  of  gravity  (CG)  to  the 
nose  tip  clockwise  around  to  vertical,  deg 


Cl 


GAMSTP 

ICHOOZ(K) 


IOPSHP 

IREDUC 

NEL 

NEMOV2 


Problem  ends  if  obliquity  angle  y  exceeds  GAMSTP, 
deg 

Selects  the  kth  plot  to  be  output  according  to  this 
table: 


ICHOOZ(K) 

PLOT 

1 

CG  axial  acceleration  versus  time 

2 

CG  axial  deceleration  versus  time 

3 

CG  lateral  acceleration  versus  time 

4 

CG  pitch  moment  versus  time 

5 

CG  velocity  versus  time 

6 

CG  Y-displacement  versus  time 

7 

CG  Z-displacement  versus  time 

8 

CG  axial  deceleration  versus  CG 
path  length,  (must  be  chosen  last 
if  chosen  at  all) 

K  =  1,2, ....NUMPLTS 

[if  NUMPLTS  -  0,  no  plots  are  output  except  two 
that  are  automatically  generated.] 

Nose  shape  selection  variable  (Figure  2.2) 

«  0  arbitrary  shape  (input  R,  THETAD,  DELS  arrays) 

1  cone 

2  cone  with  spherical  tip 

3  ogive  with  optional  cone  tip 

4  blunted  cone 

Number  of  times  to  reduce  FRAD  and  FANG  in  order  to 
complete  a  problem  with  significant  lateral  motion 
occurring  late  in  the  event  (Appendix  B) 

Number  of  longitudinal  elements  making  up  the 
projectile;  an  input  quantity  only  if  IOPSHP  -  0 

Half  the  total  number  of  circumferential  elements 
comprising  each  longitudinal  element;  must  be  an 
even  number  (Figure  2.3) 


I  •  1 


C2 


NLAY 


Number  of  target  layers  +  1  (layer  1  is  the  air 
layer  above  the  target) ;  must  be  £  5* 


NPRINT 

NSTP 

NUMNOS 

NUMPLTS 

NUMSTP 

PH1MIN 

R(I) 

RB 

RC 

RJ 

RN 

RO 


Determines  if  FREKOT  has  units  of  time  (=1)  or  path 
length  traveled  («2) 

Maximum  number  of  output  steps 

Number  of  time  steps  to  be  used  in  nose  embedment 
(see  Appendix  B) 

Number  of  plots  desired  (i.e.,  using  ICHOOZ-array  to 
select  specific  plot  types)  in  addition  to  the 
automatic  projectile  shape  and  trajectory  plots; 

0  <_  NUMPLTS  <  8 

Total  number  of  runs  to  be  generated  with  current 
set  of  data;  subsequent  runs  use  the  same  input 
except  that  DGAMAD  and  DALFAD  are  used  to  increment 
GAMMAD  and  ALPHAD,  respectively 

Angle  of  approach  4>m-tn  at  which  separation  occurs 
(i.e.,  wake  separation  angle),  deg  (see  Reference  3, 
Section  2. A) 

Average  radius  of  the  ith  longitudinal  element,  in 

Radius  of  the  conical  frustrum  (IOPSHP  »  4... blunted 
cone)  or  radius  of  conical  tip  (IOPSHP  -  3... ogive 
nose  with  optional  conical  tip),  in  (see  Figure  2.2 
for  projectile  shapes) 

Length  of  spherical  nose  (IOPSHP  *  2,  see  Figure 
2.2),  in 

Ogive  radius  (IOPSHP  -  3);  if  RJ  >_  10~12  ,  SN  is 
computed;  if  RJ  <  10-^  ,  SN  is  used  to  compute 
RJ,  in 

Sphere  radius  (IOPSHP  =2),  in 

Projectile  radius  r0  ;  used  in  time-step  computa¬ 
tions  (Appendix  B) ,  in  geometry  computations  if 
IOPSHP  t  o,  and  in  plot  heading,  in 


*  If  a  problem  requires  more  than  four  material  layers  (i.e.,  NLAY  > 
5),  then  a  parameter  variable  (NMAX)  in  PENC02D  needs  to  be  increased 
to  the  desired  NLAY-value;  the  statements  requiring  change  are  the 
first  statements  in  the  main  program  and  subroutines  SEPARATE, 

FORCES,  CALF,  BACKWARD,  and  READIN  (PENC02D  listing  in  Appendix  G) . 


SNUM(J) 


THETAD(I) 


Projectile  nose  length;  necessary  for  time-step 
computation  (Appendix  B)  and  in  geometry  computa¬ 
tions  if  IOPSHP  +  0,  in 

Young's  S-number  (Reference  2)  for  the  jth  layer; 
used  for  soil  targets;  set  =  0  for  hard  targets, 
dimensionless 

Angle  between  the  projectile  axis  and  the  line 
tangent  to  the  surface  of  the  ith  longitudinal 
element,  deg  (Figure  2.3) 


THETFD 


THETND 


TIMEF 


Projectile  aftbody  flare  angle  if  IOPSHP  >  0,  deg 
Cone  half  angle;  used  if  IOPSHP  *  1,  2,  or  4,  deg 
Case  ends  when  time  t  >  TIMEF,  sec 


TIMEI 


TITLE 


Case  begins  at  time  t  -  TIMEI,  sec 

45-character  (or  less)  title  describing  the  current 
problem;  printed  out  at  the  top  of  plots  and  printed 
output 


Impact  velocity,  in/s 

Case  ends  when  CG  velocity  becomes  less  than  or 
equal  to  VELF,  in/s 

Initial  transverse  rotation  rate  about  CG,  rad/s 


WEIGHT* 


Projectile  weight,  lb 


XICG* 


Location  of  projectile  CG  measured  from  the  nose 
tip,  in 

Location  of  aft  end  of  flare  measured  from  CG; 
positive  if  forward  of  CG,  in 

Location  of  forward  end  of  flare  measured  from  CG 
(see  XFA  above),  in 

Projectile's  transverse  mass  moment  of  inertia, 
lb-in2 


*  Appendix  F  shows  how  to  use  a  utility  program  called  MOMENT  to 

calculate  CG  location  and  moments  of  inertia  for  axisymmetric  projec¬ 
tiles  assuming  either  a  material  density  or  weight  and  an  interior 
distribution  of  material(s). 
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4 


XIRD(J) 

XLP 

YIELD(J) 

YSTOP 

ZM(J) 

ZSHIFT 


Free-surface  parameter  IT  (see  Appendix  A)  for  jth 
layer,  dimensionless 

Projectile  length,  in 

Yield  strength  for  jth  layer;  used  for  hard  targets; 
set  -  0  for  soil  targets,  lb/in^ 

Case  ends  if  |y|  (absolute  value  of  horizontal 
coordinate  of  CG)  >  YSTOP,  in 

Distance  from  ground  surface  to  the  bottom  of  the 
jth  layer  (>0),  in 

Problem  will  begin  with  projectile  nose  tip  located 
at  absolute  coordinates  Y  *  0  and  Z  =  -ZSHIFT 
(i.e.,  positive  ZSHIFT  embeds  nose  tip  downward  into 
target),  in 


ZSTOP 


Case  ends  if  Z  (vertical  coordinate  of  CG)  <  ZSTOP 


appendix  d 


LATERAL  CUTOFF 


Stable  projectiles  may  begin  with  significant  obliquity  and/or 
attack  angle,  but  will  eventually  align  themselves  with  the  velocity 
vector.  PENC02D  is  set  up  to  cut  off  the  lateral  motion  (i.e.,  set 
a  *  0,  y  =  0)  whenever  this  motion  becomes  insignificant.  For  this  to 
occur,  the  following  four  conditions  must  exist  simultaneously: 

1.  | ct {  £  0.01  deg 

2.  |y|  £  0.05  s_1 

3.  Z  £  Zcut  where  Zcut  is  a  depth  calculated  to  be  out  of  range 
of  free-surface  effects  (calculated  internally  in  the  code) 

F 

£  1.5  g's  where  =  total  lateral  acceleration 

When  these  conditions  do  occur,  the  motion  becomes  pure  axial  motion, 

indicated  in  the  printout  by  PITCH  RATE  (=y) ,  PITCH  FORCE,  and  PITCH 

_A 

MOMENT  all  equalling  zero,  with  a  very  small  (£  10  deg).  (ALPHA  is 
not  exactly  zero  because  of  the  computer  round-off  error.)  The  reasons 
for  cutting  off  the  lateral  motion  when  a  and  y  are  close  to  zero 
are  twofold: 

1.  The  projectile  begins  to  oscillate  between  small  positive  and 
negative  values  of  a  ,  with  the  frequency  increasing  in  time. 
These  oscillations  about  a  *  0  do  little  to  contribute  to  the 
overall  motion  of  the  projectile  but  do  a  lot  to  disrupt  the 
integration  scheme. 

2.  Pure  axial  motion  is  much  faster  to  integrate  because  a  number 
of  calculations  can  be  skipped  once  lateral  motion  is  cut  off 
completely. 


w 


AP PEND IX  E 

USE  OF  ELAS  TO  SMOOTH  NORMAL  STRESS  CURVE 


As  | a |  approaches  0,  the  local  normal  velocity  component  vr  for 
the  differential-area  elements  (dA's)  on  the  windward  side  of  the 
projectile  also  approaches  0  for  elements  whose  surface  is  parallel  to 
the  projectile  axis  (i.e.,  aftbody  elements  on  nonflared  projectiles). 
There  is  still  a  significant  tangential  component,  however,  as  long  as 
the  projectile  still  translates  with  some  velocity.  Thus,  the  angle 

5  -  sin-1  ^  (E.l) 


approaches  0  (i.e.,  vr  -*■  0  but  v  does  not).  This  creates  a  numeri¬ 
cal  problem  with  the  normal  stress  o  ,  because  the  o  versus  £ 
curve  is  steep  near  £  =  0.  In  fact,  the  slope  becomes  vertical  as 
£  -*■  0  . 

To  help  alleviate  the  problem,  a  smoothing  parameter  ELAS  is  intro¬ 
duced.  ELAS  defines  a  value  of  £  below  which  smoothing  will  take 
place.  The  region  0  £  <_  ELAS  of  the  a  -  £  curve  is  replaced  with 

a  curve  defined  by 


i 


r  ■. 

i  : 


o*(C)  -  o«)  •  «•» 

where  o(£)  is  the  original  a-function  for  normal  stress  prior  to 
smoothing.  Figure  E.l  shows  a  typical  a  -  £  graph  in  the  region  near 
£  *  0  for  a  6-inch  projectile  penetrating  a  soil  target  (S  *=  5)  with 
v  *  0.5Vq  ,  Vq  ■  impact  velocity  ■  500  ft/s.  The  solid  line  shows  the 
unaltered  curve.  The  dashed  curve  is  the  PENC02D  version  por  the  region 
0  <  £  <  ELAS,  demonstrating  a  smoother,  gentler  gradient  than  the 
original  curve.  The  obvious  problem  with  this  procedure  is  the  reduc¬ 
tion  in  lateral  force  experienced  by  those  differential  elements  with 
small  values  of  vn  .  Thus,  it  is  desirable  to  choose  ELAS  small  enough 

El 


-  '  "If 


'  • 


so  chat  the  total  lateral  force  Is  not  appreciably  reduced.  The  minimum 
value  of  ELAS  depends  upon  the  gradient  of  the  a  -  £  curve  and  the 
current  time  step.  If  the  value  of  ELAS  is  too  small,  one  time  incre¬ 
ment  may  be  enough  to  allow  £  to  jump  over  the  smoothing  section,  thus 
allowing  rapid  oscillations  (spikes)  to  occur  in  the  lateral  accelera¬ 
tion-time  history.  The  time-step  parameters  FRAD  and  FANG,  as  well  as 
ELAS,  must  be  chosen  so  as  to  control  spiking  without  significantly 
reducing  the  peak  lateral  acceleration.  If  ELAS  is  too  big,  the  lateral 
acceleration  will  be  attenuated.  If  both  spiking  and  lateral  peak 
reduction  occur,  a  smaller  time  step  is  needed.  Ideally,  one  should  use 
as  large  an  ELAS  value  as  possible  without  causing  attenuation. 

To  find  this  maximum  allowable  ELAS,  call  it  ELAS  ,  two  or  more 

max 

runs  should  be  executed  with  different  ELAS  values  so  that  the  lateral 
acceleration-time  curves  can  be  compared.  Figure  E.2  relates  three  runs 
made  for  the  problem  examined  in  Figure  E.l.  The  solid  curve  shows  a 
numerical  spike  at  point  C;  thus,  ELAS  used  for  that  run  must  be  too 
small.  The  alternating  dash-dot  curve  shows  reduced  peak  amplitudes  at 
points  A  and  B;  ELAS  used  for  that  run  is  therefore  too  large  (i.e., 

ELAS  >  ELAS  ).  After  a  few  runs,  an  acceptable  ELAS  ***  ELAS  can 
be  converged  upon.  In  Figure  E.2,  the  dashed  line  for  ELAS  =  0.01 
exhibits  a  very  slight  peak  reduction  with  no  numerical  spiking.  The 
conclusion  is  that  ELAS  0.01,  and  experience  with  PENC02D  has 

shown  that  ELAS  =  0.01  is  a  good  starting  value  for  most  materials. 
Nevertheless,  several  values  of  ELAS  should  be  tried  for  any  new  problem 
using  PENC02D. 


>  • 
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Normal-stress  smoothing  achieved  using  ELAS  for  6-inch  projectile  obliquely 
impacting  soil;  curve  shown  for  v  -  0.5Vo  and  z  ■  11  feet. 


of  ELAS;  oblique  impact  into  soil. 


APPENDIX  F 


MOMENT:  UTILITY  PROGRAM  TO  CALCULATE  WEIGHT,  CENTER-OF-GRAVITY 

LOCATION,  AND  MOMENTS  OF  INERTIA  FOR  AXISYMMETRIC  PROJECTILES 

When  preparing  input  for  a  PENC02D  run,  it  may  be  necessary  to 
calculate  some  or  all  of  the  projectile  parameters.  This  appendix 
presents  a  utility  computer  program  called  MOMENT  that  calculates  the 
weight,  center-of-gravity  (CG)  location,  and  polar  and  transverse  mass 
moments  of  inertia  for  any  axisymmetric  penetrator  for  which  the  weight 
(or  density),  shape,  and  location  of  each  interior  section  are  known. 

The  better  the  detail  of  the  interior  layout,  the  better  the  accuracy  of 
MOMENT'S  resulting  output  will  be.  Two  examples  will  be  set  up  to 
demonstrate  use  of  program  MOMENT. 

A  listing  of  MOMENT  is  given  in  Figure  F.l.  This  program  is  set  up 
to  run  on  a  Honeywell  DPS-1  time-sharing  system.  The  first  line  in  the 
listing  locates  additional  files  to  be  compiled  with  MOMENT  at  the  time 
of  execution;  these  files  contain  the  plot  subroutines  PLOTS,  PLOT, 
SCALIT,  SYMBOL,  AXIS13,  and  NUMBER.  All  are  routines  used  by  subroutine 
PLOTIT  to  scale  the  data  automatically,  then  draw  the  resulting  pene¬ 
trator  with  labeled  axes  on  a  Tektronix  4662  plotter  connected  to  a 
time-sharing  terminal.  There  is  a  glossary  of  input  quantities  at  the 
beginning  of  the  listing  explaining  input  variables  and  the  appropriate 
units. 

Program  MOMENT  requires  one  input  file  to  be  set  up  prior  to  execu¬ 
tion,  containing  weights  (or  densities)  and  geometry  information. 

Figure  F.2  lists  the  general  form  of  this  file.  Generally,  a  projectile 
is  divided  into  one  or  more  "figures"  for  which  the  weight  or  mass 
density  is  known.  Each  figure  must  then  be  divided  into  right-circular 
conical  frustrum  "sections"  as  shown  by  the  shaded  area  in  the  inset  in 
Figure  F.2.  For  each  figure,  the  number  of  conical  sections  (NSECT) 
must  be  specified  as  well  as  a  number  (IDENS)  that  indicates  whether  the 
weight  (IDENS  =  1)  or  density  (IDENS  =  2)  is  to  be  input.  The  next  line 
contains  the  weight  (or  density)  of  the  figure.  The  following  NSECT 
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lines  contain  the  dimensions  of  each  section  In  that  figure.  This 
Information  Is  repeated  for  each  figure  until  the  entire  projectile  has 
been  described. 

The  first  example  Is  a  blunted  conical-nosed  projectile  (Figure 

F.3)  made  of  one  material  (D6AC  steel)  with  a  known  mass  density  of 
3 

15.23  slugs/ft  .  Only  one  three-section  figure  is  needed  to  describe 
the  entire  projectile  exactly,  and  Figure  F.4  shows  the  resulting  input 
file  for  this  problem.  Line  001  indicates  the  total  number  of  figures 
to  be  used  (=1).  Line  002  shows  the  number  of  sections  ("3)  and  indi¬ 
cates  that  the  density  of  the  figure  will  be  input  (IDENS  -  2).  Line 

3 

003  is  the  mass  density  (“15.23  slugs/ft  ).  Lines  004  through  006  give 
the  required  geometry  of  all  three  sections  as  outlined  in  Figure  F.2. 
Figure  F.5  gives  the  DPS-1  time-sharing  question-answer  sequence  includ¬ 
ing  the  desired  geometric  quantities  printed  out.  (User  responses  in 
Figure  F.5  are  underlined.)  Figure  F.6  shows  the  code-generated  plot 
verifying  that  the  calculated  results  are  indeed  for  the  projectile  of 
interest  beyond  all  reasonable  doubt.  The  plot  shows  the  sectioning  of 
the  projectile  with  the  little  number  l's  indicating  the  figure  number 
to  which  the  section(s)  belong.  The  large  ")fc"  shows  the  resulting  CG 
location. 

In  July  1974,  Sandia  Laboratories  conducted  a  series  of  deep  pene¬ 
tration  tests  on  the  Watching  Hill  Blast  Range  near  Ralston,  Alberta, 
Canada.  These  experiments,  which  were  sponsored  by  the  Defense  Nuclear 
Agency,  are  documented  by  Patterson  (Reference  7).  One  of  the  penetra- 
tors  used  at  Watching  Hill  is  schematized  in  Figure  F. 7.  Because 
weights  of  individual  "compartments''  are  known  rather  than  densities, 
the  projectile  can  be  divided  into  the  eight  figures  indicated  in  Figure 
F. 7,  and  the  IDENS  ■  1  option  used.  The  resulting  input  file  is  shown 
in  Figure  F.8  and  is  self-explanatory  with  the  help  of  the  first  half  of 
sheet  1  of  Figure  F.l,  to  explain  the  terminology,  and  Figure  F.2  to 
show  the  layout  of  the  input  file.  Figure  F.9  lists  the  time-sharing 
printout,  with  user  responses  underlined,  that  results  from  executing 
MOMENT.  Finally,  Figure  F.10  plots  the  resulting  form  and  gives  the 
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OOlOXiRUN 

0020C 

0030C 

0040C 

0050C 

0060C 

0070C 

0080C 

0090C 

0100C 

OllOC 

0120C 

0130C 

0140C 

0I50C 

0160C 

0170C 

0180C 

0190C 

0200C 

0210C 

0220C 

0230C 

0240C 

0250C 

0260C 

0270C 

0280C 

0290C 

0300 

0310 

0320 

0330 

0340 

0350 

0360 

0370 

0380 

0390 

0400 

0410 

0420 

0430 

0440 

0450 

0460  75 

0470 

0480 

0490 

0500 


X;ROSDS466/CPLOTS, E; SCALE, E 

THIS  PROGRAM  CALCULATES  CG  LOCATION  AND  MOMENT  OF  INERTIA  FOR 
ANY  AXISYMMETRIC  PROJECTILE  FOR  UHICH  DENSITIES  OR  HEIGHTS  OF 
EACH  INTERNAL  COMPONENT  IS  KNOWN.  THE  PERTINENT  VARIABLES 
ARE: 


NFIGS . NUMBER  OF  (DIFFERENT-DENSITY)  FIGURES 

(I. E. SECTIONS)  WITHIN  THE  PROJECTILE. 

NSECT( I ) * I=1*NFIG£. . NUMBER  OF  SUBSECTIONS  INTO  WHICH  EACH 

FIGURE  IS  DIVIDED.  EACH  SUBSECTION 
IS  A  RIGHT-CIRCULAR  CONICAL  FRUSTRUM 
OUTSIDE  AND  INSIDE  (I.E.  INSIDE  CAN 
BE  A  HOLLOW  CONICAL  SECTION). 

IDENS( I ) , 1=1  * NFIGS . =1  TO  INPUT  FIGURE'S  WEIGHT 

2  TO  INPUT  FIGURE'S  MASS  DENSITY 

W(I) . WEIGHT  OF  I'TH  FIGURE  (LBS). 

DENSITY(I) . MASS  DENSITY  OF  I'TH  FIGURE  (SLUGS/ 

CU.FT.) 

FOR  EACH  SUBSECTION: 

R2B( J) . OUTSIDE  FORWARD  RADIUS  OF  J'TH  SUBSECTION. 

R2S( J ) . INSIDE  FORWARD  RADIUS  OF  J'TH  SUBSECTION. 

R1B(J) . OUTSIDE  AFT  RADIUS  OF  J'TH  SUBSECTION. 

R1S(J) . INSIDE  AFT  RADIUS  OF  J'TH  SUBSECTION. 

Y1(J) . DISTANCE  FROM  NOSE  TIP  TO  FORWARD  END 

OF  J’TH  SUBSECTION. 

DX(J) . THICKNESS  OF  J'TH  SUBSECTION. 

NOTE: J=1,NSECT(I)... I'TH  FIGURE. ..ALL  UNITS  ABOVE  IN  INCHES 
UNLESS  OTHERWISE  STATED. 

COPIM6N/Ai/R2B(50)iR2S(56)iRiB(50)iRis(56)lYi(56)i . 

tY2(50 ) , NTOTSEC ,ZCG, NSECT! 10 ) *DX(50) 

DIMENSION  TVOL(10),IDENS(10)*W(10)*DENSITY(10)* 

8XMASS( 10 ) > 

8E(50),F(50),G(50),H(50)*VOL(50) 

8,A(50).B(50).C(50),D(50)*X1(50)*X2(50)*MTP(50),MTT(50) 

LOGICAL  LDENS(IO) 

CHARACTER  PRTX1 
REAL  MY.MZ.MTP.MTT 
CHARACTER  X8XINPUT(2),PRTX1 
DATA  XINPUT(2)/1H;/ 

PRINT. "INPUT  FILENAME" 

READ.XINPUT(l) 

PI=3. 1416;  GRAV=32 . 2 
CALL  ATTACH! 1 >XINPUT *3*0*.) 

READ! 1,75)  LINE, NFIGS 
FORMAT! V) 

NREAD=0 

DO  10  1=1, NFIGS 
TVOL(I)=0. 

READ! 1,75)  LINE,NSECT(I).IDENS(I) 

Figure  F.l  Listing  of  program  MOMENT  (Sheet  1  of  5). 
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0510  IFdDENSd)  .E9.2)  GO  TO  11 

0520  READ(1.75)  LINE, MCI) 

0530  LDENSd>  =  .  FALSE. 

0540  GOTO  12 

0550  11  READ(1,75)  LINE.DENSITYd ) 

0560  LDENSd)  =  .TRUE. 

0570  12  IF(I.EQ.l)  GO  TO  13 

0500  NREAD=NREAD+NSECT  d-1 ) 

0590  13  CONTINUE 

0600  DO  14  J=NREAD+1 , NREAD+NSECT (I ) 

0610  READ(1,75)  LINE,R2B< J) ,R2S( J ) ,R1B( J ) ,R1S( J ) , Yl( J) ,DX( J ) 

0620  Y2(J)=Y1(J>+DX(J) 

0630  E(J)=R1B(J)-R2B(J> 

0640  F(J)=R2B(J)KY2<J)-R1B(J)XY1<J) 

0650  G( J)SR1S( J)~R2S( J) 

0660  H(J)=R2S(J)XY2(J)-R1S<J)MY1<J) 

0670  VOLU>  =  <PI/DXU>xx2>x((Y2CJ>xx3-Yl(J)xx3>X<E(J>xx2-G(J)xx2) 

0600  B/3.+(Y2(J)XX2-y1<J)kn2)k(E(J>kf<J)-G(J)kK(J))+DX(J)m 

0690  0(F( J)xx2-H( J)XX2) )/d2.**3) 

0700  14  TVOLd)-TVOLd)+VOL(  J) 

0710  10  CONTINUE 

0720  NTOTSECsJ 

0730  DO  20  I=1,NFIGS 

0740  IF(LDENSd))  GO  TO  20 

0750  XMASSd ) =U( I 7/GRAV 

0760  DENSITYd  )  = XMASSd  7/TVOL  (I ) 

0770  20  CONTINUE 

07B0  TMASS=0.;  TZCG=0.;  MN=0 

0700  DO  30  J:1,NFIGS 

0G00  IF(J.EQ.l)  GO  TO  23 

0G10  MN=MN+NSECT(J-1) 

0420  23  N=N3ECT( J) 

0030  DO  25  IC=1+MN,N*MN 

0040  XMULTsPIMDENSITY(J)/(DX(K)MX2) 

0050  DZCG=((Y2(K)kx4-Y1(K)xx4)M(E(K)MX2-G(K)xx2)/4.  ♦ 

0060  0(Y2(K)mm3-Y1(K)xx3)m(2xE(K)«F<K)>2xG(K)xh(K))/3.  ♦ 

0070  0(Y2(K)kx2-Y1(K)xx2)k(F(K)kx2-H(K)xx2)/2. )/(12.kx3) 

OOOO  DZCG=DZCGXXMULT 

0090  DMASS=DENSITY(J)XV0l(K) 

0900  TZCGSTZCG  +  DZCG 

0910  25  TMASS=TMASS  ♦  DMASS 

0920  30  CONTINUE 

0930  ZCG=TZCG/TMASS 

0940  TUEIGHT=  TMASS  X  GRAV 

0950  PRINT, 

0960  PRINT, "DO  YOU  UANT  TO  CALCULATE  MOMENTS" 

0970  READ.PRT 

0900  IF(PRT.EO.INN)  GO  TO  999 

0990  DO  40  Isl,NFIOS 

1000  DENPlsPIXOENSITY(l)/(12.XX3) 
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1018  IF(I.EQ.l)  GO  TO  43 

1021  NREAD=NREAD+NSECT<I-1> 

1030  GO  TO  44 

1040  43  NREAD=0;  MZ=0.;  MY=0. 

1050  44  CONTINUE 

1060  DO  45  J=NREAD+1 ,NREAD+NSECT( I) 

1070  X2(J)=ZCG  -  Y1CJ) 

1080  X1(J)=ZCG  -  Y2(J) 

1090  A( J)=-E< J) 

1100  B(J)=R1B<J)XX2<J)  -  R2B( J)XXI(J) 

1110  CCJ)=-G(J) 

1120  D(J)=R1S<J)XX2<J)  -  R2S( J)XXI(J) 

1130  ONE=R2B(J)KX4  +  R2B(J)X*3XR1B(J)  +  (R1B(J)XR2B(J))XM2 

1140  t  *  R2B(J)XR1B(J)XX3  +  R1B(J)XX4 

1150  TU0=R2S(J)XK4  +  R2SC J)X*3*R1S( J)  +  (R2S< J)XR1S(J))XM2 

1160  t  +  R2S(J)*R1S(J)X#3  +  RlS(J)x*4 

1170  MTPCJ)=DENPIKDX(J)x(ONE-TWO)/10. 

1180  FIVE=(X2<J)xx5-Xl(J)xx5)x(A(J>xx2-CCJ)XX2>/5. 

11*0  SIX=<X2<J)xx4-Xl<J)XX4)x(A<J)XBCJ>-CCJ>XDCJ))/2. 

1200  SEVEN=<X2<J)xx3-Xl(J)xx3)x<B(J)xx2-D(J)xx2)/3. 

1210  ADD=DENPIx<FIVE+SIX+SEVEN)/<DX<J)XX2) 

1220  MTTC J )=. 5*MTPC J )+ADD 

1230  MY=MY+MTT(J) 

1240  45  MZ=MZ+MTP ( J ) 

1250  40  CONTINUE 

1260  PRINT,;  PRINT, 

1270  PRINT  400, MY/144., MYKGRAV,MYM1. 355668/144. 

1280  PRINT, 

1290  PRINT  401, MZ/144. ,MZXGRAV,MZ«1. 355668/144. 

1300  999  PRINT, 

1310  PRINT  402, TWEIGHT.TWEIGHTX. 4536 

1320  PRINT, 

1330  PRINT  403.ZCG.ZCGX2.54 

1340  PRINT. "DO  YOU  WANT  TO  PLOT" 

1350  READ,PRT 

1360  IF(PRT.EG.IHY)  CALL  PLOTIT 

1370  400  FORMAT (2X, "TOTAL  TRANSVERSE  MOMENT  OF  INERTIA  s",lPE12.4, 

1380  *1X."SLUG-FEETxx2",/,37X,"=",E12.4,1X,"LB.-INCHE5XX2",/, 

1390  *37X,"=",E12.4,1X,"KG.-METERSXX2") 

1400  401  FORMATC2X, "TOTAL", 2X, "POLAR", 2X, "MOMENT", 2X, "OF  INERTIA", 

1410  *2X,"=",1PE12.4,1X,"SLUG-FEETXX2",/,37X,"=", 

1420  *E12.4,1X,"LB.-INCHESXX2",/,37X,"=",E12.4,1X,"KG.-METERSXX2") 

1430  402  F0RMATC2X, "TOTAL  PROJECTILE  WEIGHT-  =", 1PE12.4, 1X,"LB5.",/, 

1440  *26X,"=".E12.4,1X,"KG."> 

1450  403  F0RMAT(2X, "CENTER  OF  GRAVITY  LOCATED", 0PF10. 4, 2X, 

1460  ^"INCHES  FROM  NOSE  TIP",/,27X,"(",F9.4,2X,"CM. )") 

1470  STOP 

1480  END 
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SUBROUTINE  PLOTIT 

COMMON/A 1/R2B ( 50 ), R2S ( 50 ), RIB ( 50 ),R1S( 50), 
tYl(50) , Y2(50) ,NTOTSEC,ZCG,NSECT(1Q),DX(50) 

DIMENSION  SX(4) 

CALL  PLOTS("A") 

CALL  PL0TC2.,4.25,-3) 

SX(1)=0. 

SX(2)=0. 

RMAX=0 . 

DO  1  J=1,NT0TSEC 
SX(2)=MAX(SX(2) ,Y2(J)) 

RMAX=MAX ( RMAX , R 1 B ( J ) ) 

CONTINUE 

CALL  SCALIT(SX,7. ,2,1) 

SFX=SX<4) 

SX(2)=RMAX 

CALL  SCALIT(SX,3. ,2,1) 

SFY=SX(4) 

SFACT  =MAX< SFX, SFY) 

IF(IFIX(IFIX(SFACT)/SFACT) )  4,4,5 

ND-1 

GO  TO  7 

ND=-1 

CONTINUE 

SFINV=1. /SFACT 

H=0 .25 

CALL  SYMBOL (ZCGNSFINV.O. ,H,11,0 . »-l) 

FIRSTY=-3*SFX 

CALL  AXIS1 3 ( 0.  ,"*3. , "INCHES",*^,  .14,7 .  ,ND, 0,0 .  ,SFACT,  1 .  ,1,0) 
CALL  AXIS13( 0 . , -3 . , "INCHES" , 4 , . 14 , 4 . , ND, 1 , FIRSTY, SFACT , 1 . , 1 , 0 ) 
NF=1 

N=NSECT(1) 

DO  3  J=l,NTOTSEC 

CALL  PL0T(Y1( J)*SFINV»-R2S( J)*SFINV,3) 

CALL  PL0T(Y1( J)XSFINV,-R2B( J)»SFINV,2) 

CALL  PL0TCY2C J)XSFINV,~R1B< J)*SFINV,2) 

CALL  PL0T<Y2(J)*SFINV,-R1S(J)*SFINV,2) 

CALL  PL0T<Y1(J)*SFINV,-R2S(J)*SFINV,2) 

CALL  PL0T<Y1(J)*SFINV,R2S(J)*SFINV,3) 

CALL  PLOTCYK J)#SFINV»R2B(J)*SFINV,2) 

CALL  PL0T(Y2CJ)XSF1NV,R1BCJ)*SFINV,2) 

CALL  PL0T(Y2<J)*SFINV,R1SCJ)*SFINV,2) 

CALL  PLOT( Y1 ( J ) XSFINV , R2S( J ) *SFINV , 2 ) 

H=0.07 

H02=0.5»H 

X0=0.5»(Y1(J)+Y2(J)) 

Y0=0.25*CR1B(J)+R2B<J)+R1$<J>+R2SCJ>) 

Y00=-0 . 25X(R1B( J)+R2B( J)+R1S(J)+R2S( J)) 

IF(J.LE.N)  GO  TO  2 
NF=NF*1 


v  '-  '  .l 


r  T-  TiT. 

.■  C- 


■y.f. 
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N=N*NSECT(NF) 

2  FIG=NF 

IF(DX( J) MSFINV.LT. H)  GO  TO  3 

IF( . 5*SFINV*(R1B( J)+R2B( J)-R1S( J)-R2S( J) ) . LT.H)  GO  TO  3 
CALL  NUMBER<X0*SFINV-HO2,Y0*SFINV-HO2,H,FIG,0.,-l> 

CALL  NUMBER<X0*SFINV-HO2,Y00*$FINV-HO2, H. FIG, 0 . ,-l ) 

3  CONTINUE 

CALL  PL0TC0.,0.,999> 

RETURN 

END 
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NSECT(l) , IDENS(l) 

W(l)  or  DENSITY(l)  {depending  on  IDENS(l)} 

R2B(1) ,R2S(1) ,R1B(1) ,R1S(1) ,Y1(1) ,DX(1) 

R2B (NSECT (1) ) , R2S (NSECT (1) DX(NSECT (1) ) 
NSECT(2),IDENS(2) 

W(2)  or  DENSITY(2) 

R2B (1+NSECT (1) ) , . . .  ,DX(1+NSECT(1)) 

• 

R2B(NSECT(2)+NSECT(1)), ...  ,DX(NSECT(2)+WSECT(1)) 

« 

NSECT (NFIGS) , IDENS (NFIGS) 

W(NFIGS)  or  DENSITY (NFIGS) 

NFIGS-1  NFIGS-1 

R2B(1  +  Z  NSECT(i)),...  ,DX(1  +  Z  NSECT (i)) 

i-1  i-1 

NFIGS  NFIGS  NFIGS 

R2B(  Z  NSECT (i) ) ,  R2S(  Z  NSECT(i)) . DX(  Z  NSECT(i) ) 

i-1  i-1  i-1 


Figure  F.2  Format  of  MmENT's  input  file 
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001  1 
002  3f 2 
003  15.23 
004  I»0f3f0r0r5 
005  3»0i3>0»5i7 
006  3.1.3.1.12.8 


Figure  F.4  Input  file  for  blunted  conical-nosed  projectile. 


*0LP  MOMENT 
*FRN 

INPUT  FILENAME 
=BLNTM0M 

DO  YOU  WANT  TO  CALCULATE  MOMENTS? 
-YES 


TOTAL  TRANSVERSE  MOMENT  OF  INERTIA 


8.1801E-01  SLUG-FEET**2 
3.7930E  03  LB.-INCHESM2 
1.1090E  00  K6.-METERSW2 


TOTAL  POLAR  MOMENT  OF  INERTIA 


1.2768E-01  SLUG-FEETM2 
5.9201E  02  LB.-INCHESM2 
1.7309E-01  KB.-METERSM2 


TOTAL  PROJECTILE  HEIGHT  =  1.3255E  02  LBS. 

=  6.0124E  01  KG. 


CENTER  OF  GRAVITY  LOCATED  10.9064  INCHES  FROM  NOSE  TIP 

(  27.8540  CM.) 

DO  YOU  WANT  TO  PLOT? 

*YE8 


-  -  - 


Figure  F.5  Tine-sharing  input-output  sequence  for  blunted 
conical-nosed  projectile. 


?rrr.  ■ 

JV./.-v'.y..''. 

Ly.V.-*. 


Fll 


003  367. 

004  0*0* .47*0*0* .44 
005  .47*0*1. 04*0*. 44*1. 92 
006  1.04*0*1.6*0*2.36*2.14 
007  1.60*0*2.16*0*4.5*2.58 
008  2.16*0*2.54*0*7.08*1.76 
009  2.54*0*2.85*0*8.84*2.6 
010  2.85*0*3.02.0*11.44*1.72 
011  3. 02*0*3.16*0*13. 16*1.76 
012  3.16*0*3.22*0*14.92*1.32 
013  3. 22 *0*3. 25 *0*16. 24*9.74 
014  3. 25. .95*3. 25*1. 44*25. 98*. 2 
015  3.25*1.44*3.25*1.84*26.18* .5 
016  3. 25 >1.84* 3.25* 1.97 *26. 68* .52 
017  3. 25. 1.97*3. 25. 2. *27. 2*1. 3 
018  3. 25>2>3.25>2. >28.5>31.5 
019  4*1 
020  3. 

021  .95*0*1. 44*0*25*98*. 2 
022  1.44*0*1. 84*0*26. 18*. 5 
023  1.84*0*1. 97*0*26. 68*. 52 
024  1.97*0*2*0*27.2*1.3 
025  1*1 
026  6. 

027  2*1. 48*2. >1.48*28.5*11. 09 
028  1*1 
029  9.5 

030  2*0*2*0*39.59*7.8 
031  4*1 
032  8.5 

033  2*0*2*0*47.39*11.45 
034  1 .63*0* 1 .63*0*58.84  * .4 
035  1.63*. 85*1. 63*. 97*59. 24*. 4 
036  1.63*. 97*1. 63*. 97*59. 64*4. 6 
037  2*1 
038  3. 

039  2*1.63*2*1.63*58.84*1.16 
040  3. 25*1. 63*3. 25*1. 63*60. *.5 
041  2*1 
042  3. 

043  .85*0. *.97*0. *59. 24*. 4 
044  .97*0*. 97*0. *59.64*4.6 
045  3*1 
046  117. 

047  4. 75*3. 25*4. 75*3. 25*60. *.5 
048  4.75*1.63*4.75*1.63*60.5*3.74 
049  4. 75*0. *4. 75*0. *64. 24*1. 76 


Figure  F.8  Input  file  for  Watching  Hill  penetrator 


HOLD  MOMENT 
«FRN 

INPUT  FILENAME 
sSABOT 

DO  YOU  WANT  TO  CALCULATE  NOHENTS? 
*YES 


TOTAL  TRANSVERSE  MOMENT  OF  INERTIA  =  4.1062E  01  SLU6-FEETW2 

-  1.9040E  05  LB.-INCHES**2 

=  5.5467E  01  KG.-NETERS**2 

TOTAL  POLAR  MOMENT  OF  INERTIA  =  8.0501E-01  SLUG-FEET*#2 

=  3.7327E  03  LB.-INCHES«2 

*  1 .0913E  00  K6.-METERSM2 

TOTAL  PROJECTILE  WEIGHT  *  5.1700E  02  LBS* 

»  2.3451E  02  KG* 

CENTER  OF  GRAVITY  LOCATED  39.1027  INCHES  FROM  NOSE  TIP 

<  99.3208  CM.) 

DO  YOU  WANT  TO  PLOT? 

*YES 


Figure  F.9  Time-sharing  input-output  sequence 
for  Watching  Hill  penetrator. 


PENC02D 

WES’  TWO-DIMENSIONAL  PENETRATION  CODE 

CPENC02D 

PARAMETER  NMAX=7 

DIMENSION  RR0(70 ) ,SS0(70) ,RR1(70 ) , SSI ( 70 ) , QDECEL ( 200 ) , QSPEEDC200) 
DIMENSION  JQSET(20 ) , VAVE(3) , VPR(3) , 19(50), 9PATH(200) ,9GMP(200> 
COMMON/TRANS/SNO ,XCG0/BACK/ ALPHA/ INCR/DGAMAD, DAL  FAD, NUMSTP 
DIMENSION  ZDUM(4) ,YDUM(4) ,QM(280) ,QN(280) ,91(200) ,92(200) 
DIMENSION  TIME(200),9LTACC(200),QAXACC(200),9CGMOM(200) 
COMMON/TPLOT/TITLE,NLAY,NUMPLTS, ICH00Z(8) ,ZCUT/WT2/XIINV,XM 
COMMON/ L ABEL S/YL L ,YLR, XL ,XR,NCX, NCYL 
CHARACTER  *18YLL,XL,YLR*19,XR*19,TITLE*45 
COMMON/PXTRA/STEST , PICON , ARATIO , NEL/X2/DELPREV 
COMMON/PLTPR/FIRSTX»DELTAX, FIRSTY.DELTAY 
COMMON/BIGS/ADA(70 ) , CAT( 70 ) , SAT ( 7 0 )/FIXT/SB, SC, SE, LOGOUT 
COMMON/MSUN/ZSTOP , YSTOP , GAMSTP , VEL ,ZSHI FT 
COMMON/HOP E/VALUE(6 ) , DER( 6 ) , THETAD( 70 ) , THETND, THETFD, 

8TIMEF, TIMEI ,GAMMAD, FRE9I >NSTP, ALSTP/DOPE/W1I , FYP,NEM0V2, ALPHAD 
COMMON/DARE/WEIGHT , XICG, AGRAV , FZP , VELF 

COMMON/SHAP/DELS(70 ) ,S(70),R(70) ,THETA(70),THETAF,RB,EJ, 
tRJ , RN , RC , RO , XFA,XFF, XLP ,XCG, SN , THETAN 
COMMON/FINK/ZM( NMAX) , GAMI , RADDEG, DEGRAD, R90 , R180 , R270 , R360 , IOPSHP 
C0MM0N/DAN4/XSALP(24),XCALP(24)/R0UND/ALAS,ELAS 
COMMON/NEWT/FRAD, FANG.DTMAX, DTMIN, FREKOT , IREDUC,NUMNOS,NPRINT 
COMMON/INTERM/AKASER( 70 , NMAX) , BKBSER(70 , NMAX) , GKCATE(NMAX) , 

AHRD1 ( NMAX ) , HRD2 ( NMAX ) 

COMMON/YOUNG/LMEDIA( NMAX) , SNUM(NMAX), SHRGD( NMAX) , DENSITY! NMAX) , 
*YIELD(NMAX)/SEPARE/SYNEFEE,TANFEE,PHIMIN,RCAV(70),S0(10) 

LOGICAL  LOGDT , LMEDIA, LPRNT , LOGOUT 
REAL  MU 

DATA  MU/1 . 633E4/ , BET/0 . 1744/ , GAM/90 . 28/ 

DATA  AK/1 -  253/ , BK/3 . 67/ , CK/5 . 9/ 

FS9(J,X)=S9RT(X(J)**2+X(J+1)**2) 

AGRAV-386 . 07 

R360=6. 2831853;  R180=0.5#R360;  R90=R360*0 .25;  R270=R360«0 . 75 
RADDEG-360 ./R360 
DEGRAD=R36 0/360. 

SZE8=5. 

CALL  PLOTS! , ,3) 

CALL  PLOT(0 . ,8. ,-3) 

LMEDIA(1)=.T. 

SNUM(1)=1.E10 

ZM(1)=0. 

1000  URITE(6,4) 

CALL  READIN 
FRADSV=FRAD 
FANGSV=FANG 
VCHKSV- ,2KVEL 

NUMB-0; *SN0=SN;  XCG0=XCG;  IBACK=1 

AKA-AKMMU 

BRB=BK*S9RT(BET) 

GKC=CKNQAM 
PICON-R180/NEMOV2 
IF(ELAS)  15,15,16 

15  ELAS:0 .01 

16  ALA$sR90/ELAS 
IMULT=NEM0V2/2 

1001  GAMI SGAMMADN  DEGRAD 
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GULF*(GAMMAD-ALPHAD)*DEGRAD 
THETAN=THETND*DEGRAD 
THETA F=THETFD* DEGRAD 
IF(IOPSHP.NE.O)  GO  TO  103 
XNL=NEL 

DO  107  1=1, NEL 
THETAC I ) =THETAD( I ) *DEGRAD 
107  CONTINUE 

103  S( 1 )=XCG  -  DELS(l)/2. 

NUMB=NUMB  +  1 

I F ( NUMB . GT . NUMSTP )  GO  TO  1000 

XCG=XCGO 

SN=SN0 

FRAD=FRADSV 

ICNT=0 

VCHEK-VCHKSV 

FANG=FANGSV 

YITO=XCG*SINCGAMI) 

ZITO=-XCG*COSCGAMI)  -  ZSHIFT 
LOGCUT=. F. 

DYITO=-VEL*SINCGULF) 

DZITO=  VELXCOS(GULF) 

IFCNUMB . EQ . 1 )  GO  TO  104 
IF( IBACK . EG . 1 )  GO  TO  S35 
CALL  BACKWARD 
IBACK=1 
GO  TO  B35 

104  IF(IOPSHP.NE.O)  GO  TO  106 
DO  105  1=2, NEL 

IM1=I-1 

S(I)=S(IMl)-( (DELS( IM1 )+DELS( I ) )/2 . ) 

105  CONTINUE 
GO  TO  200 
106  CONTINUE 

STEST  =  XCG  -  XLP 
DELS(2)=DELS(1) 

CALL  CONSET 

GOTO(110,120,130,140), IOPSHP 
110  CALL  SHPKS1000 ) 

GO  TO  201 
120  CALL  SHP2 
GO  TO  201 
140  CALL  SHP4 
GO  TO  201 

130  CALL  SHP3(01000) 

201  INL=NEL 

NEL=NEL  ♦  1 

S(NEL)=S(INL  )-DELS(  IND/2.-.0005 
DELS(NEL)=. 001 

R(NEL)=(R(INL)  *  DELS( INL )«TAN(THETA( INL ) )/2. )/2. 

THETA (NEL )=-ATAN(R( NEL)/. 0005) 

200  CONTINUE 

CALL  SEPARATE 
WRITE(6,2)  IOPSHP, NEL 

2  FORMATCl  SHAPE  IS",I2,4X,"NEL=*»,I3,//,4X,"I", 
09X,"R(I)",9X,"S(I)",7X,"THETA(I)",6X,"DELS<I)",6X,"RCAV(I)") 


G3 


PENC02D 

WES'  TWO-DIMENSIONAL  PENETRATION  CODE 

41  CONTINUE 

RMAXs0 . ;  NLM=4*NEL+1 ;  NEX=NLM-1 

D0210I=1,NEL 

J2=2#I 

JI=J2  -  1 

K1=NLM  -  J2 

K2=NLM  -  J1 

D0VE2=DELS(I)/2. 

QM(J1)=S(I)  +  D0VE2 
QM( J2)=S( I )  -  D0VE2 
QM(K2)=QM( J1 ) 

QM(K1)=QM( J2) 

DR=D0VE2*TAN( THETAC I )  ) 

QN( J1 )=R( I )  -  DR 
QN<J2)=R(I>  +  DR 
QN(K2)S**QN(  J1 ) 

QN(K1)=-QN(J2) 

RMAX= AMAX1 ( RMAX , R ( I ) ) 

X=THETA(I)MRADDEG 

WRITE<6,5>  I,R(I),S(I),X,DELS(I>,RCAV(I) 

55  CONTINUE 
RR0(I)=R(I) 

SS0(I)=SCI) 

J=NEL-I*1 

RRKI)sR(J) 

210  SSlCI)s_S(J) 

VALP=QM(1>  -  QM(2*NEL) 

VBET=2 . XRMAX 

IF(VALP. LT. VBET)G0TO212 

CALL  SCALIT(QM,SZE8*NEX» 1) 

QN(NEX+1 ) =QM( NEX+1 ) 

QN(NEX+2)=QM(NEX+2) 

G0T0213 

212  CALL  SCALIT(QN,SZE8.NEX,1) 

QM( NEX+1 ) =QN( NEX+1 ) 

QM(NEX+2)=QN(NEX+2) 

213  CONTINUE 
NCX=25 

XL=18H  DIST  FROM  CG  < 

XR=3HIN) ;  NCYL=21 
YLL-18HX-SECTI0N  RADIUS  < 

YLR=3HIN) 

SZD-SZE8+ . 15 

CALL  SYMBOLd . ,  SZD, . 2 , 16HPR0JECTILE  SHAPE, 0., 16) 
CALL  PLOTIT(NEX,9N,QM,0,SZE8> 

9N(NEX+1)=0. 

DL=PICON 
DO  10  J=1,IMULT 
ALP  =  (J-.5)*DL 
XSALPCJ)  =  SINCALP) 

10  XCALPCJ)  =  COS( ALP) 

DO  11  J=IMULT+1 'NEM0V2 
JJ=NEM0V2+1-J 
XSALP(J)  =  XSALP(JJ) 

11  XCALP(J)  =-XCALP(JJ) 

DO  915  JJ=1,NLAY 
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IF( .NOT.LMEDIA(JJ))  GO  TO  918 
DO  916  1=1, NEL 

AKASER(I,JJ)=AKA/(2.XR<I)*SNUM(JJ>> 

916  BKBSER<I,JJ)=BKB*S9RT(2.*R(I)>/SNUM<JJ> 

GKCATEC JJ) =GKC/SNUM( J J ) XX2 
GO  TO  915 

918  HRD1(JJ)=6.3298»YIELD(JJ) 

HRD2( J J ) =4 . 74735XSQRT ( DENSITY( J J ) *YI ELD( J J ) ) 

915  CONTINUE 

DO  829  1=1, NEL 
ATEH  =  R90  -  THETA(I) 

CAT ( I >  =COS ( ATEN ) 

SAT(I)=SIN(ATEH) 

829  ADA(I)=DLkR(I)xDELS(I)/COS(THETA(I)) 

835  WRITE(6,4> 

VALUE(l)=YITO 

VALUE(2)=ZIT0 

VALUE(3)=DYITO 

VALUE(4)=DZIT0 

VALUE(5)=GAMI 

VALUE(6)=M1I 

VB=VEL 

W=UEIGHT ;  I9SET=1 
JQSET(1)=1 
JQSET(2)=NSTP*1 
CALL  CALF 
ABALF=ABS( ALPHA) 

PLAST=0. 

MXX=0 

IF( SN . L E . 0 . )  SN=0.5*XLP 
S  EGL  EN =SN/NUMNOS 
IBACK=1 

NCNT1=0;  TIND=TIMEI;  LOGDT=.F. 

NF=lJ  FNF=FREKOT 
LPRNT=.F. 

IFCNPRINT .EG. 2)  LPRNT=.T. 

DEL=SEGLEN/VB 
PATH=0.;  NBODY=l 

ptst=nbodyxfreqi*xlp 

PTDIF=PTST|  IQ<1)=1 
JG=1 

CALL  FORCES(VALUE, TIMED 

WRITE(6,3)  JQ,YITO,ZITO, GAMMAD, NCNT1 , DYITO, DZITO ,1411 , 
8MXX, T IMEI , FYP , FZP , SB , DEL , FYP/U , FZP/W, SB/U, ALPHA , SC , 

SSC/U  SE  W 

91(1 )=YITO;  92(1)=ZIT0;  QGMP(1)=GAMI-R180;  TIME(1)=TIMEI 
9LTACC(1)=SB/W;  9AXACC(1)=SC/U;  9CGM0M(1)=SE 
9DECEL(1)=-9AXACC(1);  9SPEED(1)=VB/12. I  9PATH(1)=0. 
9LTMIN=9LTACC( 1 ) ;  9LTMAX=9LTACC< 1 ) ;  9DECMX=9DECEL(1) 

DO  400  J9=2,NSTP 
220  CONTINUE 

IF(LOGDT)  GO  TO  227 
IF(PATH.GT.SN)  GO  TO  502 
DT1=SEGLEN/VB 
GO  TO  504 

502  DT1=FRADKR0/VB 


PENC02D 

WES'  TWO-DIMENSIONAL  PENETRATION  CODE 

504  ABV6=ABS(VALUE(6) ) 

IF( ABV6 . GT .0.5)  GO  TO  501 

DEL=DT1 

GO  TO  503 

501  DT2=FANG/ABV6 
DEL=MINCDT1,DT2) 

503  DEL=MINCDEL,DTMAX> 

DEL -MAX (DEL > DTMIN) 

227  CALL  FORCES( VALUE*  TIND) 

VB=FSQ(3* VALUE) 

IF(VB.LE.VELF)  MXX=7 
IFCTIND..GE.TIMEF)  MXX=3 
CALL  CALF 
ABALF=ABS( ALPHA) 

IFCABALF.GE. ALSTP)  MXX=9 

IFCABALF.LT. 90.)  GO  TO  224 

CALL  BACKWARD 

IBACK=-IBACK 

IQSET=IQSET+1 

JQSET(IQSET)=JQ 

GO  TO  223 

224  IF(MXX.GT .2)  UO  TO  223 
IF(LPRNT)  GO  TO  20 
IF(TIND.LT.FNF)  GO  TO  45 
GO  TO  22 

20  IFCPATH.LT. FNF)  GO  TO  45 

22  NF=NF+1 

FNF=NFMFREKOT 
GO  TO  223 

45  IFCLOGCUT )  GO  TO  288 

IFCVALUEC2) . GT.ZCUT)  GO  TO  288 

IF(ABS(VALUE(6)).GT. 0.05. OR. ABALF.GT. 0.01)  GO  TO  288 

IFCABSCSB/U) .GT. 1 .5)  GO  TO  288 

LOGCUT=.T. 

VALUE(3)=-VB*SIN( VALUEC5) ) 

VALUEC4 ) =  VBMCOSCVALUEC5)) 

VALUEC6)=0. 

GO  TO  227 

288  IFCLPRNT)  GO  TO  24 

IFCTIND+DEL .GT . FNF)  DEL=FNF-TIND 
24  VPRCD-VALUEC3) 

VPR(2)SVALUE(4) 

VPR(3)=VALUE(6) 

186  DERC3)=DEL»FYP*XM 
DERC4)=DEL*FZPXXM 
IFCLOGDT)  GO  TO  187 

IF(VB.GT.VELF*2.0XFSQ(3»DER))  GO  TO  187 
LOGDT-.T. 

DEL=DELPREV 
DELSAV-DEL 
GO  TO  186 

187  IFCLOGCUT)  GO  TO  48 
DERC6 )-DELKSEKXIINV 
GO  TO  47 
48  DER(6)-0. 

47  DO  239  11=1*3 
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II2sII+2 

IP=II 

IFCII.NE.3)  GO  TO  238 
112=6 
IP=5 

238  VALUEC II2)=VALUECII2)  +  DERCII2) 

VAV£CII)=CVALUECII2)  *  VPR(II))/2. 

DERCIP)=VAVECII)*DEL 

239  VALUEL IP) = VALUEC IP)  +  DERCIP) 

TIND=TIND  +  DEL 

NCNT1=NCNT1+1 

DELPREV=DEL 

PATH=PATH  ♦  FSQC1,DER) 

IFCPATH.LT. PTST)  GO  TO  220 
NBODY=NBODY  ♦  1 
PLAST=PATH 

PTST=NBODY#FREQI«XLP 
IQCNBODY)=JQ 
GO  TO  220 

223  IFCJQ.EQ.NSTP)  MXX=2 

IFCABSCVALUEC1)) .GT.YSTOP)  MXX=4 
IFC VALUE C2). LT .2STOP )  MXX=5 
I FC VALUEC 2) . GT . MAXCXCG, XLP-XCG) )  MXX=6 
V5= VALUEC 5 )*RADDEG 
IFCV5.GT. GAMSTP)  MXX=8 

WRITEC6 ,3)  JQ, VALUEC 1 ) , VALUEC 2 ), V5 , NCNT1 , VALUEC 3 ), VALUEC 4 ) , 
8VALUEC6 ) ,MXX, TIND. FYP , FZP , SB , DELPREV , FYP/W, FZP/W, 
tSB/W, ALPHA . SC , SC/W , SE , U 

3  FORMAT  C/» IX, 6HSTEP  =  .  14 ,  UX,  13HY-P0SITI0N  =,E13.6, 

84X, 12HZ-P0SITI0N  =  , E13 . 6 , 4X, 14HPITCH  ANGLE  =,E13.6, 

84X,  9HCOUNT  =, I4,/,22X, 13HY-VELOCITY  =,E13.6,4X, 

0 12HZ-VEL0CITY  =,E13.6,4X,14HPITCH  RATE  =,E13.6,4X, 
89HRETURN  =,I4./,1X,6HTIME  =, El 1 . 4 , 4X, 13HY-F0RCE  =  , 

AE13.6,4X, 12HZ- FORCE  =, E13 .6 ,4X, 14HPITCH  FORCE  =,E13.6, 

84X,  7HDEL  =,E11.4,/,22X,5HY-G'S»7X, 1H=»E13.6»4X»5HZ-G'S» 
06X, 1H-, E13.6,4X, 9HPITCH  G' S,4X, 1H= , EI3. 6 , 4X, 7HALPHA  =, 

8E11 .4*/»22X» 13HAXIAL  FORCE  =  , E13 . 6 , 4X , 12HAXIAL  G’S  =  , 
tE13.6,4X, 14HPITCH  MOMENT  =  ,1PEU .4,2X, 

812HPROJ  WEIGHT=, Ell .4) 

IFCLOGCUT)  GO  TO  19 
IFCVB.GT . VCHEK)  GO  TO  19 
IFCICNT.GE. IREDUC)  GO  TO  19 
VCHEK=0 . 5* VCHEK 
FRAD=FRAD/3. 

FANG=FANG/3. 

ICNT=ICNT+1 
19  CONTINUE 

IFCLOGDT)  DEL=DELSAV 
QIC JQ)=VALUEC1) 

Q2CJQ)=VALUEC2) 

QGMPC JQ) = VALUEC 5) -R180 
TIMEC JQ)=TIND 
QLTACCC JQ)=SB/U 
QAXACCC  JQ)=SC/kJ 
QCGMOMC JQ ) =SE 
QDECELC JQ)=-QAXACCC JQ) 


r 


PENC020 

WES'  TWO-DIMENSIONAL  PENETRATION  CODE 


QSPEED<JQ)=VB/12. 

QPATHC JQ)=PATH 
NCNT1=0 

QLTMAX=MAX(QLTMAX,QLTACC( JQ) ) 

QLTMIN=MINCQLTMIN,QLTACC(JQ)) 

QDECMX=MAX( QDECMX , QDECEL ( JQ ) ) 

IFCMXX.GT.2)  GO  TO  450 
400  CONTINUE 
450  NBODY=NBODY  +  1 
IQ(NBODY)=JQ 

IF(NBODY. EQ . 2 . OR . PATH-PLAST . GT . 0 . 55*PTDIF)  GO  TO  255 

NBODY=NBODY  -  1 

IQ(NBODY)=JQ 

255  JQSETCIQSET+l )=JQ+3 

IF( JQ.LT.4)  GO  TO  777 

WRITE<6,420)  QPATH< JQ ) , QPATHC JQ)/12 . , CTIMEC JQ)-TIMEI >*1000 . , 
*QSPEED( JQ) >QLTMIN, QLTMAX, QDECMX 

420  F0RMAT(//3X, "FINAL  STATISTICS", /3X,16(1H*), /IX, "PATH  LENGTH  ", 
*"TRAVELED=",E13.5,"  IN.  (",E13.5,"  FT. )",/lX,"TIME  ELAPSED=", 
8E13.5,"  MSEC", /IX, "FINAL  SPEED=",E13.5,"  FPS",/1X,"MIN.  LATERAL  • 
*," ACCELERATION5", El 3. 5,"  G'S",/IX,"MAX.  LATERAL  ACCELERATION5", 
8E13.5,"  G'S",/1X,"MAX.  AXIAL  DECELERATION-", El 3. 5,"  G'S") 

T ADD= AMAX1 ( XLP , 2#R0 ) 

CALL  SCALIT(Q1 ,S2E8 , JQ, 1 ) 

X1=Q1<JQ+1);  DX1=Q1( JQ+2) 

CALL  SCALIT(Q2,SZE8, JQ, 1) 

X2=Q2< JQ+1) ;  DX2=Q2( JQ+2) 

63  CONTINUE 

ZDUM(1)5Q2(1);  ZDUM(2)=Q2(1);  YDUM(1)=Q1(1);YDUMC2)=Q1(1) 

DO  234  J=1,JQ 

ZDUM( 1 ) =AMIN1 ( Q2( J )-TADD, ZDUM( 1 ) ) 

ZDUM(2)=AMAX1 (Q2( J)+TADD,ZDUM(2) ) 
YDUM(l)sAMINl(QI(J)-TADD,YOUMa)) 

234  YDUM(2)=AMAX1(Q1( J)+TADD,YDUM(2) ) 

CALL  SCALIT(ZDUM,SZE8,2, 1) 

CALL  SCALITtYDUM. SZE8,2, 1 ) 

X1=YDUM(3) 

DX1=YDUM(4) 

X2=ZDUM( 3) 

DX2=ZDUMC4) 

Q2( JQ+1)  =  X2 

Q2( JQ+2)  5  DX2 

Ql< JQ+1)  =  XI 

QIC JQ+2)  =  DX1 

DELSCAL  -  AMAX1 ( DX1 , DX2 ) 

IF(DELSCAL-DXl)  350,352,350 
352  Q2( JQ+2)  5  DELSCAL 

CALL  SETUP CX2 , SZE8 , DEL SCAL , ISTRT , DX2 , DX1 ) 

Q2( JQ+1 )  5  ISTRT  #  DELSCAL 
GO  TO  354 

350  QKJQ+2)  5  DELSCAL 

CALL  SETUPIXl , SZE8 , DELSCAL , ISTRT, DXi , DX2 ) 

QKJQ+l)  =  ISTRT  «  DELSCAL 
354  YLl=18H  Z  ( 

XL  =18H  Y  ( 

YLR=3HIN);  NCYL-36 
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XR-YLR;  NCX=NCYL 
FIRSTX  =  Q1<JG+1) 

DELTAX  =  9HJ9+2) 

FIRSTY  =  Q2CJQ+1) 

DEITAY  =  92U9+2) 

DO  333  J=1 .NBODY 
JUSE=I9(J) 

INDEX=2 

330  IFUUSE.LT. JQSET(INDEX))  00  TO  33X 
INDEX=INDEX  ♦  1 

GO  TO  330 

331  INDEX=INDEX  -  1 
INDY=< INDEX/2) *2 

IF< INDY. E9 . INDEX)  60  TO  447 

CALL  PLOTPROC  91UUSE)  ,  92( JUSE) ,9GMP( JUSE)  , SSO ,RRO ) 
GO  TO  333 

447  CALL  PLOTPRO  (Q1  (JUSE)  ,Q2(  JUSE)  .QGMPUUSE) ,  SSI,  RR1) 
333  CONTINUE 

NCHS=-FIRSTY/DELTAY 

XNC=NCHS 

CALL  AXIS13(0. ,XNC,1H  , 1,0. ,SZE8,-1,0, 0 . ,0. , 1 . , 1,0) 

XST=- . 5 

IZE=SZE8 

DO  125  LZ-l.IZE 

XST=XST+1.0 

CALL  CALCMPCXST , XNC, 0 ,  1 ) 

125  CALL  CALCMP(XST,XNC-.07,16,1) 

J0=1 

123  J0=J0*1 

IFCJO.GT.NMAX)  GO  TO  124 
IF(ZM( J0)+ZDUM(1) )  , 124 , 124 
YST=ZM(JO)/DELTAY 
YST=XNC-YST 

CALL  CALCMP<0.,YST,0,1) 

CALL  CALCMP(SZE8,YST,15,1) 

GO  TO  123 

124  CALL  PLOTIT(JQ.92,Q1,0,SZE8) 

IF( NUMPLTS . E9 . 0 )  GO  TO  1000 
XL=18H  TIME  ( 

XR=4HSEC);  NCX=34 

CALL  SCALIT(TIME,SZE8, JQ, 1) 

DO  149  J*l, NUMPLTS 

GO  TO  <131, 132. 133, 134, 135, 136, 137, 138), ICHOOZ(J) 

C *****  PLOT  AXIAL  ACCELERATION  VS  TIME  ***** 

131  YLL-18H  AXIAL  ACC  < 

YLR=4HG'S);  NCYL=29 

CALL  PLOTIT( J0.9AXACC.TIME, 1 ,SZE8) 

GO  TO  149 

C*****  PLOT  AXIAL  DECELERATION  VS  TIME  ***** 

132  YLL=18H  AXIAL  DEC  < 

YLR=4HG'S);  NCYL=29 

CALL  PL0TIT(JQ,QDECEL,TIME,1,SZE8) 

GO  TO  149 

CKNMMK  PLOT  LATERAL  ACCELERATION  VS  TIME  ***** 

133  YLL*18H  LATERAL  ACC  < 

YLRs4HG’S) )  NCYL=27 
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CALL  PLOTIT! JQ.QLTACC.TIME.l *SZE8) 

GO  TO  149 

C*****  PLOT  CG  MOMENT  VS  TIME  ***** 

134  YLL=18H  CG  MOMENT  C 

YLR=6HLB-IN);  NCYL=31 

CALL  PLOTIT! JQ*QCGM0M*TIME*1*SZE8) 

GO  TO  149 

C*****  PLOT  VELOCITY  VS  TIME  ***** 

135  YLL=18H  VELOCITY  ( 

YLR=4HFPS);  NCYL=30 

CALL  PLOTIT(JQ* QSPEED*  TIME* 1  * SZE8) 

GO  TO  149 

C*****  PLOT  Y-DISPLACEMENT  VS  TIME  ***** 

136  YLL=18H  Y  ( 

YLR=3HIN> }  NCYL=36 

CALL  PLOTIT!JQ*Ql , TIME* 1  * SZE8) 

GO  TO  149 

C**X**  PLOT  Z-DISPLACEMENT  VS  TIME  ***** 

137  YLL=18H  Z  ( 

YLR=3HIN>;  NCYL=36 

CALL  PL0TIT(JQ*Q2 , TIME* 1  * SZE8 ) 

GO  TO  149 

ex****  PLOT  AXIAL  DECELERATION  VS  PATH  LENGTH  ***** 

138  YLL=18H  AXIAL  DEC  ( 

YLR=4HG'S) ;  NCYL=29 
XL=18H  PATH  LENGTH  ( 

XR=3HIN);  NCX=26 

CALL  SCALIT!QPATH*SZE8* JQ.l) 

CALL  PLOTIT! JQ*  QDECEL .QPATH* 1  * SZE8) 

149  CONTINUE 
777  CONTINUE 

IF! DAL FAD)  1003,1002,1003 

1002  IF(DGAMAD)  1003.1000.1003 

1003  GAMMAD=GAMMAD  +  DGAMAD 
ALPHAD=ALPHAD  ♦  DALFAD 
GO  TO  1001 

4  FORMAT! 1H1) 

5  FORMAT! 2X, I3,5X,F10.5,4F13.5) 

END 


Cxxxx*  CALCULATE  CAVITY-RADIUS  ARRAY  AND  SEPARATION  POINT!S)  ***** 
SUBROUTINE  SEPARATE 

COMMON/SEPARE/SYNEFEE*  T ANFEE, PH,RCAV! 70 ) *  SO! 10 ) 

8/PXTRA/DUM! 3) , NEL/SHAP/DELS! 70 ) * 

&S!70)*R!70)*TH!70), THETAF*  DUMN! 5) * RO*XFA,XFF*XLP*XCG,SN 
RCAV! 1>=R! 1 ) 

J  =2 
1=0 

10  IF(PH)  15»15»11 

11  IF(PH.GT.TH( J) )  GO  TO  20 

15  RCAV! J)=MAX!R! J) *RCAV! J-1)) 

J=J  +  1 

IF! J.GT.NEL)  GO  TO  100 
GO  TO  10 
20  1=1+1 
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S0(I)=S(J)+0.5*DELS( J) 

R0=R(J)-0.5*DELS(J)*TAN(TH(J)) 

R0SQ=R0**2 

RTF=2*R0*TANFEE 

22  RCAV(J)=SQRT(R0SQ+RTF*(S0(I)-S(J))) 
IF(RCAV(J) . LT . R(  J ) )  GO  TO  15 
J=J  +  1 

IF(J.GT.NEL)  GO  TO  100 
GO  TO  22 
100  1=1+1 

SO ( I )=-l . E20 

RETURN 

END 


C****x  DETERMINE  AXES  STARTING  LOCATION  FOR  TRAJECTORY  PLOT  xxxxx 
SUBROUTINE  SETUP (X, SIZE, DSCALE, ISTRT , D1 , D2) 

ISTRT  =  X/SIZE 
ISTRT  =  ISTRT  -  4 
15  IF(ISTRTXDSCALE-X)  20,20,10 
20  ISTRT  =  ISTRT  +  1 
GO  TO  15 

10  ISTRT  =  ISTRT  -  1 

IF(D1.LT. (.75*02))  ISTRT  =  ISTRT  -  1 

RETURN 

END 


Cxxxxx  PLOT  PROJECTILE  SHAPE  AT  (YCG.ZCG)  ON  TRAJ.  PLOT  xxxxx 
SUBROUTINE  PLOTPRO(YCG,ZCG,GAMMAP,S,R) 

COMMON/PLTPR/'FIRSTX,  DELTAX, FIRSTY, DELTAY 
COMMON/PXTR A/DUMMY ( 3 ) , NEL 
DIMENSION  S(1),R(1),X(150),Y(150) 

SGMP=SIN( GAMMAP ) 

CGMP =COS ( GAMMAP ) 

DO  10  J=1 , NEL 
ONE=S( J ) *SGMP 
TWO=S( J ) XCGMP 
THREE=R(J)XSGMP 
FOUR=R( J )*CGMP 
X(J)=ONE  +  FOUR  +  YCG 
Y( J)  =  ZCG  -  TWO  +  THREE 
JJ  =  2*NEL  -  J  +  1 
y(JJ)=ONE  -  FOUR  +  YCG 
10  \(JJ)=ZCG  -  (TWO  +  THREE) 

NPTS  =  2*NEL  +  1 
X(NPTS)  =  X(l) 

Y(NPTS)  =  Y( 1 ) 

X(NPTS+1 )  =  FIRSTX 
X(NPTS+2)  =  DELTAX 
Y(NPTS+1 )  =  FIRSTY 
Y(NPTS+2 )  =  DELTAY 
CALL  LINE(X,Y,NPTS,1,0,12) 

CALL  SYMBOL ( ( YCG-FIRSTX)/DELTAX, (ZCG-FIRSTY)ZDELTAY, .07,11 


Gil 
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*,0. ,-l) 
RETURN 
END 


Cxxaxx  DRAW  And  LABEL  AXIS  kxkxx 

SUBROUTINE  AXIS13(X,Y, IBCD,NC,H, SIZE,NN, IXY,XMIN,DX, 
C  SPACE, ITIC.NP10) 

COMMON/P AXIS/ I BEGIN 
DIMENSION  IBCD(l) 

DATA  IBCD1/1H  / 

NCHARS=ABS(NC) 

XY=IXY;  ANGLE=XYX90. ;  TIC=ITIC 
SPAC=SPACE;  G=H;  XT=X;  YT=Y 
GVERT=G;  LINE=2 
SIZ=SIZE;  P10=NP10 
IF(SIZ)  1,16,100 

1  LINE=3 
SIZ=-SIZ 

100  ND=NN 
NA=NN 

IF(ND)  3,2,2 

2  NDIG=ND  +  1 
GO  TO  30 

3  ND=0 
NDIG=0 

30  NSPACE=SIZ/SPACE  +  .5 
FNSPAC=NSPACE 
TL=FNSPACXSPAC 

IF(YT  +  TLXXY  -  29.)  4,4,31 

31  CALL  SYMBOL (XT , YT+ . 1 , . 14 , 13HAXIS  TOO  LONG, 90. ,13) 
RETURN 

4  H02=.5XG 

H07=H02/3.5 

TEMP=-NP10 

POWER=10.x*TEMP 

DEL X= POWER *DX;  XNUMB=POWERXXMIN 

ANG=1.  -  XY 

ALAB=1. 

IF(NC)  5,6,6 

5  ALAB=-1. 

6  TICXY=(2.XXY-1. )XALABXTICX.1 
XTIC=TICXYXXY;  YTIC=TICXYxANG 
TICM1 = . 16  -  TICX.05 

LIH=3;  NMAX=0;  J=0 
DO  12  I=J,NSPACE 
NDIGIT=NDIG 
CALL  PLOT(XT,YT,LIN) 

LIN=LINE 

CALL  PLOT(XT+XTIC,YT+YTIC,2) 

IF(G)  7,11,7 

7  TEMP=-ND 

ITEMP=ALOG10(ABS(XNUMB)+.5X10.XXTEMP)  +  1. 

IF(ITEMP)  70,70,71 
70  ITEMP=1 
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71  NDIGIT=NDIGIT+ITEMP 

IF( ABS(XNUMB) .LT. 1 .  E-6)  XNUMB=0. 

IF(XNUMB)  9,10,10 

9  NDIGIT=NDIGIT+1 

10  IF(NDIGIT-NMAX)  77,77,76 

76  NMAX=NDIGIT 

77  FDIGIT=NDIGIT 
CENTER=C7.XFDIGIT-3. JXH07X.5 
XAN0=-ALABXXY«( TICM1+CENTER)  -  CENTER 
YANO=ALABKANG*(TICM1+H02>  -  H02 

CALL  NUMBERCXT+XANO,YT+YANO,G,XNUMB,0.,NA> 
XNUMB=XNUMB+DELX 

11  CALL  PLOT (XT » YT,3) 

XT=XT+SPACXANG 

12  YT=YT+SPACXXY 
IF(NCHARS.GT.l)  GO  TO  20 
IF( IBCD( 1 ) . EQ . IBCD1 )  RETURN 

20  ANOLOC=IBEGINXG*ALAB 

ANOCHR-( SIZ  -  NCHARSXGVERTJX.5 

XST=X  -  XY*ANOLOC  +  ANGXANOCHR 

YST=Y  +  ANGXANOLOC  ♦  XYXANOCHR 

CALL  SYMBOL (XST,YST,G, IBCD, ANGLE, NCHARS) 

IF(NPIO.EQ.O)  RETURN 

ANOLOC-( IBEGIN  +  ALABX . 5)XGXALAB 

ANOCHR= ( SIZ  *  GVERTX(40-NCHARS))X.5 

XST=X  -  XYXANOLOC  ♦  ANGXANOCHR 

YST=Y  +  ANGXANOLOC  +  XYXANOCHR 

CALL  NUMBERLXST , YST , 5 . XH07 , P10 , ANGLE, -1 ) 

16  RETURN 
END 


Cxxxxx  CONNECT  PIECES  OF  AXIS  LABEL  TOGETHER  xxxxx 
SUBROUTINE  TLABtKP, ICHOSE, TITLE, NC) 

COMMON/L ABELS/YLL , YLR,XL ,XR,NCX, NCYL 
CHARACTER  X18YLL ,XL , YLRX19.XRX19 
CHARACTER  LLABxlB , RLABXI9, TITLEX42 
CHARACTER  FMT1X17 , FMT2X20 , FMT3X17 

DATA  FMT1/17HC A18,"10  ",A19)/,FF1T2/20H(A18,"10  ",A19,"  ">/, 

A  FMT3/17H( A18, A19,"  ")/ 

IFCICHOSE. EQ.2)  GO  TO  11 
LLAB=XL;  RLAB=XR;  NC=NCX 
GO  TO  13 

11  LLAB=YLL;  RLAB=YLR;  NC=NCYL 
13  IFCKP)  21,22,23 

22  NSP=0 

ENCODEt TITLE, FMT3)  LLAB,RLAB 
GO  TO  25 

23  IFCKP-9)  24,24,21 

24  NSP=4 

ENCODE ( TITL E, FMT2)  LLAB.RLAB 
GO  TO  25 
21  NSP=5 

ENCODECTITLE, FMT1 )  LLAB.RLAB 

25  NC=NC+NSP 
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RETURN 

END 


C*»*x*  DETERMINE  AXIS  ANNOTATION  STARTING  LOCATION  ***** 

SUBROUTINE  ANOSTRT  C LOCATE, KLOG, Q1 , QSIZE, DELQ,RINC,NN) 
POWER=10.**C-KLOG) 

ANOMIN=Ql*POUER 

ANOMAX=ANOMIN  +  QSIZE*DELQ*POUER/RINC 
AMINAB=ABS( ANOMIN) 

AMAXAB-ABSCANOMAX) 

MXLOG=0 ;  MNLOG=0 

IF( AMINAB .GT .(1.0))  MNLOG=ALOG10(AMINAB> 

IF( AMAXAB.GT. ( I . 0 ) )  MXLOG=ALOG10( AMAXAB) 

L0CATE=5  +  MAXCMXLOG, MNLOG)  +  NN 
IF(MXLOG-MNLOG)  10,10,20 
10  IF( ANOMIN . LT . ( 0 . ) )  LOCATE=LOCATE  +  1 
GO  TO  30 

20  IFCANOMAX.LT. CO.))  LOCATE=LOCATE  +  1 
30  RETURN 
END 


C*****  PLOT  YARRAY  VS  XARRAY  ***** 

SUBROUTINE  PLOTITC JQ, YARRAY, XARRAY, ISCALE.SZ) 

COMMON/P AXIS/ I BEGIN 
DIMENSION  XARRAYC 1 ) , YARRAYC 1 ) 

CHARACTER  TITLE*42 , AAA*S/1H  / 

IFCISCALE.EG.O)  GO  TO  30 
CALL  SCALITCYARRAY,SZ, JQ, 1 ) 

30  X1=YARRAYC JQ+1) 

XDEL  =YARRAY  C  JQ+2 ) 

KP=0 

IFCXDEL . GT . 0 . )  KP=AL0G10CXDEL ) 

IFCKP.LT. 0)  KP=KP-1 
CALL  TLABCKP,2,TITLE,NC) 

CALL  ANOSTRTCLOCATE, KP ,X1 , SZ,XDEL , 1 . , 1 ) 

IBEGIN=LOCATE 

CALL  AXIS13C0. , 0 ., TITLE, NC, . 0 9, SZ, 1 , 1 ,X1 , XDEL , 1 . , 1,KP) 
XSTRT =- . 5 
DO  10  L=1,SZ 
XSTRT=XSTRT+1. 

CALL  CALCMPC 0 . ,XSTRT, 0 , 1 ) 

10  CALL  CALCMPC. 07, XSTRT, 16,1) 

X1=XARRAYC  JQ+1 ) 

XDEL=XARRAYC JQ+2) 

KP=0 

IFCXDEL. GT. 0. )  KP=ALOG10 CXDEL ) 

IFCKP.LT. 0)  KP=KP-1 
CALL  TLABCKP,1, TITLE, NC) 

IBEGIN=5 

CALL  AXIS13C0. ,0. , TITLE, -NC, . 09, SZ, 1 , 0 ,X1 ,XDEL , 1 . , 1 , KP) 
XSTRT=- . 5 
DO  20  L-l , SZ 
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XSTRT=XSTRT+1. 

CALL  CALCMPCXSTRT  ,  0  . ,  0 , 1 ) 

20  CALL  CALCMP(XSTRT, .07,16,1) 

CALL  LINE!XARRAY,YARRAY,JQ,1,0,12> 
EMCODEC AAA , 55 )  TITLE 
55  F0RMATCA8) 

IF(AAA. EQ.8H  DIST)  GO  TO  25 
IF(AAA . EQ .8H  DIST  F)  GO  TO  25 
CALL  HEADING 

25  CALL  CALCMP!0.,0. ,0000,2) 

RETURN 

END 


CXXXXX  DETERMINE  SCALE  FOR  ARRAY  GIVEN  SZ  AXIS  LENGTH  xxxxx 
SUBROUTINE  SCALIT!ARRAY,SZ,NPTS,INC> 

DIMENSION  ARRAY(1),TIS(9),VM(9) 

DATA  TIS/. 07918, .17609, .30103, .39794, .47712, .60206, 
8.69897, .77815, .90309/, VM/1. 2, 1.5, 2. 0,2. 5, 3. 0,4.0, 
85.0,6.0,8.0/ 

TMX=0.  ;TMN=0. 

DO  10  1=1 , NPTS, INC 
TMX= AMAX1 C  TMX , ARRAY(I)) 

TMN=AMIN1 ( TMN , ARRAY( I ) ) 

10  CONTINUE 
SZUSE=SZ 

50  SUSE=(TMX-TMN)/SZUSE 
IF  (SUSE)  12,118,12 
12  QLOG=0. 

IF  (SUSE.GT.O. )QLOG=ALOG10(SUSE) 

KLOG=QLOG 
CHAR=KLOG 
TISSA=QLOG-CHAR 
VMULT=1.0 
IF(TISSA >5,110,20 
5  TISSA=TISSA*1. 

KLOO=K»,3G-1 

I  FT  TT'wSA .  LE.  1 .  E-10  )  GO  TO  110 
20  DO  41  J=1 , 9 

IF  !TISSA.GT.TIS! J))GO  TO  41 
VMULT=VM(J> 

GO  TO  110 
41  CONTINUE 

KLOG=KLOG*l 

110  ARRAY! NPTS*2)SVMULTX 10 . KXKLOG 
XTEST=TMN/ARRAY(NPTS+2> 

NTESTsXTEST 

IFtNTEST . LE.XTEST)  GO  TO  120 

NTEST=NTEST-l 

GO  TO  120 

118  ARRAY!NPTS*2>=0. 

NTEST=0 

120  ARRAY! NPTS*1 >sNTESTX ARRAY! NPTS+2) 

RM„H=ARRAY!NPTS*1 > 

RMAX=RMIN*SZX ARRAY! NPTS+2 ) 
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I F  ( TMX . GT . RMAX . OR . TMN . L T . RMIN ) 60  TO  60 
RETURN 

60  SZUSE=SZUSExfl.95 
60  TO  50 
END 


Cxxxxx  PLOT  DESCRIPTIVE  HEADING  ON  FINAL  PLOTS  xxxxx 
SUBROUTINE  HEADING 
COMMON/ SHAP/DMY( 286) ,RO/TP LOT/TITLE 

COMMON/MSUN/DUMN ( 3  ) , VEL/DOPE/DUN<  3 ) , ALPHAD/DARE/WEIGHT 

COMMON/HOP E/DUM2 ( 86 ) , GAMMAD 

CHARACTER  BADX2/2HD=/ , BBDX8/8H  IN. ,  W=/, 

IBCDX8/8H  LBS,  V=/, BDDX4/4H  FPS/,BEDX6/6HALPHA=/, 
8BFDX13/13H  DEG.,  GAMMA=/,BGDX5/5H  DEG./,TITLEX45 
SZZ=6. 

SZF=SZZ-.2 
SZG=SZF- . 2 

CALL  SYMBOL ( . 35, SZZ, .1, TITLE, 0. ,45) 

Rl=2 . KRO 

CALL ’ SYMBOL ( .35,SZF, . 1 , BAD,  0  .  ,2) 

CALL  NUMBER ( 999. , SZF, . 1 , R1 , 0  .  ,2) 

CALL  SYMBOL ( 999 . , SZF , . 1 , BBD, 0  .  ,8) 

CALL  NUMBER< 999. , SZF, .1, WEIGHT, 0. ,2) 

CALL  SYMBOL ( 999 . , SZF , . 1 , BCD,  0  .  ,8 ) 

R1=VEL/12. 

CALL  NUMBER ( 999 . , SZF, .1,R1,0.,-1) 

CALL  SYMBOL ( 999 . , SZF, .1,BDD,0.  ,4) 

CALL  SYMBOL ( . 35 , SZG, . 1 , BED, 0 . , 6 ) 

CALL  NUMBER ( 999 . , SZG, . 1 , ALPHAD, 0 . ,2) 

CALL  SYMBOL ( 999. , SZG, .1,BFD,0. ,13) 

CALL  NUMBERC 999., SZG,. 1, GAMMAD, 0.,-l) 

CALL  SYMBOL ( 999 . , SZG, . 1 , BGD, 0 . , 5 ) 

RETURN 

END 


Cxxxxx  GEOMETRY  GENERATOR  FOR  CONICAL  NOSE  xxxxx 
SUBROUTINE  SHP1(X) 

COMMON/SHAP/DELSt  70 ) ,S(70)»R(70),THETA(70)»THETAF,RB»EJ 
*RJ,RN,RC,RO,XFA,XFF,XLP,XCG,SN,THETAN 
TTHETF=TAN(THETAF) 

IFCTHETAN-l.E-4)  1,1,2 

1  IF(SN-l.E-4)  99,99,3 
3  C3=R0/SN 

THETAN=ATAN(C3) 

GO  TO  5 

2  C3=TAN(THETAN) 

SN=RO/C3 

5  C1=XCG-SN 
IFINISH=0 
1  =  1 

10  IF(I.NE.l)  S(I)=S(I-l)-((DELSCI-l)+DELS(I))/2.) 

IF(S(I)  .LT.CDGOTOIOO 


PENC02D 

WES’  TWO-DIMENSIONAL  PENETRATION  CODE 


Rd  )  =  (XCG-Sd )  )*C3 
THETA(I)=THETAN 
GO  TO  150 

100  IF(S(I) .  LE.XFF.  AND.Sd) . GE.XFA)  GO  TO  120 
Rd  )=RO 
THETA(I)=0. 

GO  TO  150 
120  CONTINUE 

Rd  )  =  RO  +  (XFF  -  S(I))  KTTHETF 
THETA(I)=THETAF 
150  IFdFINISH.EO.l)  RETURN 
CALL  CONVERL IFINISH.I) 

GO  TO  10 
99  MRITE(6 »  98) 

98  FORMAT (//» IX, "SHAPE  1  THETAN  AND  SN  ARE  ZERO.") 
RETURN  1 
END 


Cxxxxx  GEOMETRY  GENERATOR  FOR  CONICAL  NOSE  TOPPED  BY  A  SPHERE  xxxxx 
SUBROUTINE  SHP2 

COMMON/SHAP/DELS( 70 ) »S(70)#R(70),THETA(70),THETAF»RB»EJ» 
*RJ,RN,RC,R0,XFA,XFF,XLP,XCG,SN,THETAN/X2/DELPREV 
TTHETF=TANCTHETAF) 

TTHETN=TAN( THETAN) 

C1=XCG-SN 

C4=XCG-RC 

RNSQ=RNXRN 

C=RN-RC 

RB=SQRT(RNSQ-CXC) 

IFINISH=0 

1=1 

SUMI=0. 

10  IFd.NE.l)  Sd)=Sd-l)-<<DELSd-l)+DELSd))/2.) 

SUMIM1  =  SUMI 
SUMI=SUMI+DELS( I ) 

IF(Sd) .LE.C4)GOT0120 
IF(I.GT.l)GOTO110 
XL=SQRT(RNSQ-(RN-DELSd  )  )  ) 

THETA(l)=ATAN(XL/DELSd)) 

C=RN-DELS( 1 )X . 5 
R(1)=SQRT(RNSQ-CXC) 

GO  TO  150 
110  CONTINUE 

C=RN-SUMIM1 

RFI=SQRT<RNSQ-CXC) 

C=C-DELS(I)x.5 

R(I)=SQRT(RNSQ-CKC) 

C=RH-SUMI 

RAI=SQRT(RNSG-CXC) 

THETACI )=ATAN( (RAI-RFI )/DELS(I) ) 

GO  TO  150 
120  CONTINUE 

IF(Sd).LT.Cl)GOTO130 

THETAd)=THETAN 
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RC I ) =RB+(C4-SC I ) >/TTHETN 
CO  TO  150 
130  CONTINUE 

IF(S( I ) . LE.XFF. AND.SCI) .GE.XFA)G0T0145 

R(I)=RO 

THETACI)=0. 

GO  TO  150 
145  CONTINUE 

RC I ) =RO*ABS ( S( I ) -XFF) XTTHETF 
THETA(I)=THETAF 
150  I FC I FINISH. EQ.  1)  RETURN 
II=I 

CALL  CONVERCIFINISH,I) 

IFC I . EQ . I I )  SUMI=SUMI-DELPREV 
GO  TO  10 


CXXXXX  GEOMETRY  GENERATOR  FOR  OGIVE  WITH  OPTIONAL  CONE  TIP  CRB,SN>0)  XX 
SUBROUTINE  SHP3(X) 

COMMON/SHAP/DELS(70 ) ,SC70),RC70), THETA( 70 ) » THETAF, RB , EJ , 
«RJ,RN,RC,RO,XFA,XFF,XLP,XCG,SN,THETAN 
IFCRJ.LT. l.E-12)  GO  TO  110 
PHI1=ARSINCEJ/RJ) 

T3=RJXC0S(PHI1)-R0 

PHI2=ARC0S(T3/RJ) 

SN2=RJXSIN(PHI2)  -  EJ 
IFCRB.LT. l.E-12)  SN=SN2 
GO  TO  130 

110  IFCSN.LT. l.E-12)  GO  TO  120 
SN2=SN 
RB=0. 

PHI3=ATANCRO/SN) 

T4=SNx . 5/COS CPHI3) 

T5=CEJ+SNK.5)/SINCPHI3> 

RJ=SQRTCT5XT5*T4XT4) 

GO  TO  130 
120  CONTINUE 
URITEC6,1) 

RETURN  1 

1  FORMATC1X. "SHAPE  3  RJ  AND  SN  ARE  ZERO") 

130  CONTINUE 
C1=XCG-SN 
TTHETF=TANCTHETAF) 

RJRJ  =  RJ  X  RJ 
T8  =  SN  ♦  EJ  -  XCG 
IFIN1SH=0 
1  =  1 

IFCRB.LT. l.E-12)  GO  TO  14 
PHIRB=ARCOSCCT3+RB)/RJ) 

XX=SN-RJXSINCPHIRB) 

TTHETN=RB/XX 
THETAN=ATANC  TTHETN) 

IFCDELSCD.GT.C0.8XXX))  GO  TO  11 
NC0NE=IFIXCXX/DELSC1))+1 
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DXX=XX/NCONE 
THETA(1)=THETAN 
DELS( 1 ) =OXX 
XXX=DXX/2. 

S<I)=XC6-XXX 
Rtl)=TTHETNXXXX 
12  1=1+1 

XXX=XXX+DXX 
DELS( I ) =DXX 

I F( I . GT . NCONE)  GO  TO  10 
SI I ) =XCG-XXX 
RII)=TTHETNXXXX 
THETA(I)=THETAN 
GO  TO  12 
11  DELS! 1 ) =XX 

SI  1  ),=XCG-XX/2 . 

THETA! 1)=THETAN 
RU)=RB/2. 

1=2 

10  SI I )=S! I-l)-t I  DELS! 1-1 )+DELS( I ) )/2 . ) 

14  IFISII) . LE.C1)  GO  TO  140 
T77=T8+SI I ) 

T7=SQRTI RJRJ-T77XT77 ) 

R< I )=T7-T3 
C=T77-DELSt I )X . 5 
RTH1=SQRTIRJRJ-CXC) 

C=C+DELSI I ) 

RTH2=SQRTtRJRJ-C*C) 

THETAII )=ATAN! IRTH1-RTH2)/DELS!1 ) ) 

GO  TO  150 

140  IFISII). LE.XFF. AND. SII). GE.XFA)  GO  TO  145 
RII) =RO 
THETAII )=0 . 

GO  TO  150 
145  CONTINUE 

RII) =RO+ABSISII )-XFF)XTTHETF 
THETAII )=THETAF 
150  IFIIFINISH.EQ.l)  RETURN 
CALL  CONVERI I FINISH « I ) 

GO  TO  10 
END 


Cxxxxx  GEOMETRY  GENERATOR  FOR  BLUNTED  CONICAL  NOSE  xxxxx 
SUBROUTINE  SHP4 

COMMON/SHAP/DELSI70 )»SI70)»R(70),THETAI70) , THETAF, RB , EJ , 
IRJ » RN,RC,RO,XFA,XFF,XLP,XCG»SN » THETAN 
TTHETN=TAN(THETAN) 

TTHETF=TANITHETAF) 

C  XCG  MUST  BE  GREATER  THAN  NOSE  LENGTH 
RI1)=0.5XRB 

THETA! 1)=ATANIRB/. 001) 

DELS! 1 ) =0 . 001 
St  1 ) =XCG- .00  05 
RC=  RB/TTHETN 
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SNT  =  RO/TTHETN 
SN=  SNT-RC 
Cl=  XCO-SN 
IFINISH-0 
1  =  1 

GO  TO  150 

10  SCI )=S( 1-1 )-( (DELSCI-1 )+DELS( I ) )/2 . ) 
IF(S(I).LT.C1)  GO  TO  142 
THETA(I)=  THETAN 
R(I )=(XCG-S(I ) )#TTHETN*RB 
GO  TO  150 
142  CONTINUE 

IF(S(I ) . LE.XFF. AND.S(I) . GE.XFA)  GO  TO  144 
R( I )=  RO 
THETACI )=0. 

GO  TO  150 
144  CONTINUE 

Rd)  =  RO  ♦  (XFF  -  SCI))  KTTHETF 
THETA(I)=THETAF 
150  IFCIFINISH.EQ.l)  RETURN 
CALL  CONVERC IFINISH, I ) 

GO  TO  10 
END 


C*«**«  LONGITUDINAL  ELEMENT  THICKNESS  GENERATOR  (CONSTANT  L/U  FOR  DA) 
SUBROUTINE  CONSET 

COMMON/SHAP/BEIS(70),SC70),R(70),THETA(70),THETAF,RB,EJ, 
*RJ,RH,RC,RO,XFA,XFF,XLP,XCG,SN,THETAN 
COMMON/PXTRA/STEST .PICON , ARATIO , NEL/X2/DELPREV 
LOGICAL  LOG, L0G2 

FLOC(I,X)=S(I>-0.5*DELS(I)+X-XCG 

ISET=0 

SUZE=1 . 

SGNM=1. 

L0G=.F. 

L0G2=.F. 

ISET=0 

RETURN 

ENTRY  CONVER (I FINISH, I) 

IF( LOG)  GO  TO  5 
IFCI.EQ.l)  GO  TO  2 
IFCISET.EQ.l)  GO  TO  8 
DELPREV=DELS(I) 

DEL SC  I ) =PICON#R ( I ) NCOS ( THETACI ) ) X ARATIO 
TEST=ABSC  CDELPREV-DELSC I ) )/DELPREV ) 

IFCTEST.GT. .003)  RETURN 

ISET=1 

RETURN 

8  IFCL0G2)  GO  TO  1 
SUZEP=SUZE 
SGNMP=SGNM 
SUZE=FLOC(I,SN) 

SGNM=SIGNC 1 . , SUZE) 

IFCSGNMP-SGNM-1.)  1,5,3 
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3  DELS( I)=FLOC( 1-1 ,SN) 

LOG=.T. 

RETURN 
5  LOO=.F. 

L062= . T . 

1  IF(FLOC(I,XLP>.LE.O.)  GO  TO  100 

2  1=1+1 
ISET  =  0 

IF(I.GT.70)  GO  TO  102 
DELS!I )=DELS( 1-1 ) 

RETURN 

100  SLENG=0 . 

DO  101  J=1,I-1 

101  SLENG=SLENG+DELS! J) 

DELS( I )=XLP-SLENG 
NEL=I 

IFINISH=1 

RETURN 

102  URITE(6,103) 

103  FORMAT!//, 3X, "SUBROUTINE  CONVER ... I>7 0 ...  DECREASE  NEM0V2  OR  ", 

A  "INCREASE  ARATIO .") 

STOP 

END 


Cxxxxx  CALCULATE  COMPRESSIVE  NORMAL  STRESS,  SUM  FORCES  A  MOMENTS  xxx 
SUBROUTINE  FORCES! VALUE, TIND  ) 

PARAMETER  NMAX=7 

COMMON/YOUNG/ LMEDI A (NMAX) , SNUM(NMAX) , SHRGD(NMAX) , DENSITY(NMAX) , 
AYIELD(NMAX)/INTERM/AKASER!70,NMAX) ,BKBSER(70,NMAX) ,GKCATE(NMAX) , 
AHRDl(NMAX) ,HRD2(NMAX) 

LOGICAL  LMEDIA, LOGCUT 
COMMON/PXTRA/DUMN  C  3 ) , NEL/BACK/ALPHA 
COMMON/FIXT/SB, SC, SE, LOGCUT 

COMMON/SEPARE/SYNEFEE, TANFEE»PHIMIN»RCAV(70 ) , SO ( 10 ) 
C0MM0N/DAN4/XSALP! 24) ,XC ALP (24 ) /ROUND/ ALAS, EL AS 
COMMON/ BIGS/ A DA( 70 ),CAT(70),SAT(70) 

COMMON/DOPE/W1 I , FYP , NEM0V2 , ALPHAD 

COMMON/FINK/ZM(NMAX),GAMI,RADDEG, DEGRAD, R90, R180 ,R270,R360 , IOPSHP 
COMMON/DARE/WEIGHT, XICG, AGRAV, FZP, VELF 

COMMON/SHAP/DELS( 70 ) ,S(70),R(70),THETA(70),THETAF,RB,EJ, 

ARJ , RN, RC, RO.XFA ,XFF,XLP,XCG, SN, THETAN 
DIMENSION  VALUE! 1 ) 

TINDX=TIND 

Z=VALUE(2) 

DY=VALUE( 3 ) 

DZ=VALUE(4) 

VCG=SGRT ( DYXDY+DZXDZ) 

ALF=ALPHAXDEGRAD 
SALF=SIN( ALF) 

SGAM=SIN! VALUE! 5) ) 

CGAM=COS! VALUE(S) ) 

U1B=VALUE(6) 

YB=  CGAMXDY  ♦  SGAMxDZ 
YC=-SGAMXDY  ♦  CGAMXDZ 
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TDVZ=0.5*VALUE(6)/YC 

SB=0. 

SC=0 . 

SE=0. 

EPS=VCG*SYNEFEE 

JL=0 

JLP1=1 

DO  410  IS1,NEL 
SI  =  SCI) 

RI  =  RC I ) 

1  IFCSI.6T.S0CJLP1))  GO  TO  2 
JL=Jl+l 

JLP1=JLP1*1 
GO  TO  1 

2  IF(JL.EQ.O)  GO  TO  3 
RM=TDVZ*CS0(JL)-SI)**2 
C2=<SI-XCG)*SALF 
RC2=RCAV(I)**2 
RI2*RI**2 

3  FA  =  Z  +  CGAM*SI 
Q  =  SGAMXRI 
BTEST=MIN(FA+Q, FA-Q) 

IF(BTEST)  10,50,50 

50  IF(VALUE(5)-4. 71238)  412,412,410 
10  DA=ADA(I ) 

YBI  =  YB  -  W1B*SI 
CATE  =  CAT ( I ) 

SATE  =  SAT ( I ) 

RCTHCG=RI*CGAM*CATE 
32  FF=U1B*RI 

TB=0 . ;  TC=0 . ;  TE=0. 

DO  411  J=1 , NEM0V2 
SALP  =  XSALP(J) 

CALP  =  XCALP(J) 

208  ZRN=FA-Q*CALP 

IFCZRN.GE. 0. )  GO  TO  411 
DO  12  K=2,NMAX 
IF(ZM(K)+ZRN)  12,13,13 

12  CONTINUE 
GO  TO  411 

C*****  CHECK  FOR  FREE  SURFACE  EFFECTS  ***** 

13  ZDAMAG=FA-Q*CAIP*(1.+SHRGD(K)*SATE)+SHRGD<K)*RCTHCG 
DO  14  KDAM=K-1,K+1 

IF(ZM(KDAM)+ZDAMAG)  14,15,15 

14  CONTINUE 

15  IF(K. EQ.KDAM)  GO  TO  19 
IF(LMEDIA(K) )  GO  TO  17 
IF(LMEDIACKDAM) )  GO  TO  18 
IF(YIELD(K) . GT .YIELDCKDAM) )  GO  TO  18 
GO  TO  19 

17  IF< .NOT.LMEDIA(KDAM) )  GO  TO  19 
IF(SNUM(K) . GT.SNUM(KDAM) )  GO  TO  19 

18  K=KDAM 

19  CONTINUE 

C*k**k  END  FREE  SURFACE  EFFECT  CHECKS  ***** 

B23=-CALP*SATE 
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QB:YBI 

QC=YC  -  FFXCALP 
ZC=B23XQB+CATE#9C 
IFCZC.LE.O. )  GO  TO  411 
ZCC=YCXCATE  ♦  YBXB23 
IFtZCC.LT.O.)  GO  TO  411 
Cxxxxx  CHECK  FOR  SEPARATION  xxxxx 
IF(JL.EQ.O)  GO  TO  430 
IF(ZCC.GE.EPS)  GO  TO  430 
C1=-RI*CALP  -  C2  +  RM 
C3=SQRT(RC2  -  RI2XSALPXX2) 

IF(CALP)  21,21,22 

21  IF(C1.LT.C3)  GO  TO  411 
GO  TO  430 

22  IF(-C1.LT.C3)  GO  TO  411 
Cxxxxx  END  CHECK  FOR  SEPARATION  xxxxx 

430  ZB=SALPXQB 

B21=-CALPXCATE 

ZA-B21X0B-SATEXQC 

SMALL V=S9RT(ZAXX2+ZBXX2*ZCXX2) 

VNV2=ZC/SMALLV 
VNV=SQRT( VNV2) 

IFCLMEDIA(K))  GO  TO  451 
SIGMAS=(HRD1(K)4-SMALLVXHRD2(K))XVNV 
GO  TO  453 

451  SIGMAS=(AKASER(I,K)+BKBSER<I,K)XSHALLV-ZRNKGKCATE(K>k 
«  VNVIXVNV 

453  AA=1. 

BB=ARSIN( VNV2) 

IF(BB.LT.ELAS)  AA=SIN(BBXALAS) 

Y3--AAKSIGMASXCATE 
TC=TC  ♦  Y3 

IF(LOGCUT)  GO  TO  411 
Y2=AAXSIGMASXSATEKCALP 
TB=TB  ♦  Y2 
TE=TE  -  CALPXY3 

411  CONTINUE 
TC=TCXDA 
SC=SC*TC 

IF(LOGCUT)  GO  TO  410 

TB=TBXDA 

SB=SB  +  TB 

DATE=DAXRIXTE 

SE=SE  +  (DATE  -  TBXSI) 

410  CONTINUE 

412  SC=2.XSC 
SB=2.XSB 
SE=2.XSE 

FYP=CGAMXSB-SGAMXSC 

FZP=SGAMXSB+CGAM«SC 

RETURN 

END 


Cxxxxx  CALCULATE  ATTACK  ANGLE  ALPHA  xxxxx 
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SUBROUTINE  CALF 

PARAMETER  NMAX=7, NMAX1=NMAX+1 

COMMON/HOPE/ VALUE( 6  >/BACK/ALPHA/FINK/DUMMY(NMAXl ) , RADDEG, 
8DEGRAD, R9Q  >  R180 ,R270,R360 
IF(VALUE(3) )  515,510,515 
515  IF( VALUE(4) )  525,520,525 
525  TH=ATAN<-(VALUE(3)/VALUE<4))) 

GO  TO  530 

510  IF(VALUE(4).LT.O.)  GO  TO  512 
ALPHA=VALUE(5)*RADDEG 
GO  TO  35 

512  ALPHA=( VALUE(5)-R180)*RADDEG 
GO  TO  35 

520  IFCVALUE(3).LT.O. )  GO  TO  522 
ALPHA=(VALUE(5)-R270 )*RADDEG 
GO  TO  35 

522  ALPHA=< VALUE(5)-R90 )*RADDEG 
GO  TO  35 

530  IF<VALUE(4)>  531,520,533 
533  IF(VALUE(3).LT.O. )  GO  TO  532 

ALPHA=CVALUE(5)-TH-R360)*RADDEG 
GO  TO  35 

532  ALPHA=<VALUE<5)-TH)*RADDEG 
GO  TO  35 

531  ALPHA=(VALUE(5)-TH-R180)*RADDEG 
35  IF(ALPHA-180. )  31,30,30 

31  IF( ALPHA+180 . )  32,32,34 

32  ALPHA=ALPHA  +  360. 

GO  TO  31 

30  ALPHA=ALPHA  -  360. 

GO  TO  35 
34  RETURN 
END 


CXXKXX  REVERSE  GEOMETRY  FOR  TUMBLING  PROJECTILE  CALPHA>90  DEG.)  *** 
SUBROUTINE  BACKWARD 
PARAMETER  NMAX=7 , NM3=NMAX*3 

COMMON/TRANS/SNO ,XCGO/HOPE/VALUE( 6 )/BACK/ALPHA/FINK/DUMMY(NM3 ) , 
8R90 ,R180/PXTRA/STEST , DL , ARATIO , NEL/B IGS/ ADA ( 7  0 ) ,CAT(70),SATC70) 
COMMON/SHAP/DELS( 70 ) ,S(70),R(70),THETA(70), DUM( 10 ) ,XCG, SN 
DATA  ICNT/0/ 

DIMENSION  SBCK(70),RBCK(70),DLBCK(70),THBCK(70) 

ICNT=ICNT  +  1 
ITEST=<ICNT/2)X2 

VALUE(5)=VALUE(5)  -  SIGN( 1 , VALUEt 5)-R180 ) XR180 
ALPHA=ALPHA  -  SIGNC 1 , ALPHA ) *180 
IF(ICNT.EQ.ITEST)  GO  TO  15 
SN=0 .001 

XCG=-S<NEL)  +  0.5*DELS(NEL) 

GO  TO  16 

15  SN=SN0 
XCG=XCG0 

16  DO  10  J=1 ,NEL 
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SBCK! J)=-S! JN) 

RBCK(J)=R(JN) 

DLBCK( J)=DELS( JN) 

10  THBCK! J)=-THETA!JN) 

DO  20  J=1 ,  NEL 
S!J)=SBCK! J) 

R! J)=RBCK!J) 

DELS! J)=DLBCK( J) 

THETA!J)=THBCK! J) 

ATEH=R90  -  THETA(J) 

CAT( J ) =COS( ATEH ) 

SAT! J)=SIN!ATEH) 

20  ADA!J)=DLXR!J)XDELS! J)/COS!THETA!J)) 

100  FORMAT!//, "PROBLEM  RESTARTED  —  PROJECTILE  GOING  BACKWARD", 
»"  NOW",///) 

110  FORMAT!//, "PROBLEM  RESTARTED  —  PROJECTILE  GOING  FORWARD", 
*"  NOW",///) 

CALL  SEPARATE 

IF!ICNT.EQ.ITEST)  GO  TO  35 
WRITE!6,100) 

GO  TO  36 

35  WRITE!6,110) 

36  RETURN 
END 


Cxxxxx  READ  AND  ECHO-PRINT  PROBLEM  INPUT  xxxxx 
SUBROUTINE  READIN 
PARAMETER  NMAX=7 

COMMON/TPLOT/TITLE»NLAY,NUMPLTS» ICH00Z!8) , ZCUT 
COMMON/YOUNG/ LMEDI A ! NMAX ) , SNUM!NMAX) , SHRGD!NMAX) , DENSITY! NMAX) 
t , YI ELD! NMAX) /ROUND/ ALAS, EL  AS 

COMMON/SEPARE/SYNEFEE, TANFEE, PHIMIN , RCAV! 70 ) , SO! 10 )/WT2/XIINV,XM 
LOGICAL  LMEDI A 
DIMENSION  XIRD!NMAX) 

CHARACTER  X1DDUM/1H  /.TITLEX45 

COMMON/PXTRA/STEST , PICON, ARAT 10, NEL/INCR/DGAMAD, DALFAD,NUMSTP 
COMMON/NEWT/FR AD , FANG , DTMAX , DTMIN, FREKOT , IREDUC , NUMNOS , NPRINT 
COMMON/MSUN/2STOP , YSTOP , GAMSTP, VEL  >ZSHIFT 

COMMON/FINK/ZM! NMAX), GAMI,RADDEG, DEGRAD, R90 , R180 , R270 , R360 , IOPSHP 
COMMON/SHAP/DELS! 70 ) ,S!70),R!70),THETA!70),THETAF,RB,EJ, 

BRJ , RN,RC, RO,XFA ,XFF,XLP , XCG, SN , THETAN 
COMMON/DARE/WEIGHT , XICG , AGRAV , FZP , VEL  F 

COMMON/HOPE/VALUE!6) ,DER!6) ,THETAD!70) ,THETND,THETFD,TIMEF, 
*TIMEI,GAMMAD,FREQI,NSTP,ALSTP/D0PE/W1I,FYP,NEM0V2,ALPHAD 
111  FORMAT ! 1H1 ) 

1000  FORMAT ! 1216 ) 

1001  FORMAT! 4E18 . 6 ) 

1002  FORMAT ! 5E14 . 4 ) 

998  FORMAT ! 3X , A45 ) 

READ!5,998,END=999)  TITLE 
997  FORMAT!10X,A4  >) 

WRITE! 6 , 997 )  TITLE 
WRITE!6,901) 

901  FORMAT!//, 3X,37!1HX),/,3X, 
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*"*  WES  2D  NEWTONIAN  PENETRATION  CODE 
*/,3X,37!lH*)) 

READ!5, 1000 )  NEM0V2 , IOPSHP ,NSTP»NLAY» NUMPLTS, NUMSTP , NUMNOS 
H.IREDUC.NPRINT 

WRITE! 6 , 2000 )  NEM0V2, IOPSHP,NSTP,NL AY, NUMPLTS, NUMSTP, NUMNOS 
t , IREDUC ,NPRINT 

2000  F0RMAT!/,3X,"NEM0V2=",I4,2X,"I0PSHP=",I4,2X,”NSTP=",I4,2X, 

*"NLAY=",I4,2X,"NUMPLTS=",I4,2X,«NUMSTP=",I4,2X,"NUMN0S=",I4,2X, 
i"IREDUC="  >I4,2X,"NPRINT=",I4) 

IF! NUMPLTS. EQ.  0)  GO  TO  950 
READ(5,1000)( ICHOOZC I ) , I =1 , NUMPLTS ) 

WRITEC6, 2002) (ICHOOZC I), 1=1, NUMPLTS) 

2002  FORMAT!/, 3X, "ADDITIONAL  PLOTS . . . " , 7 ( 14  )  ) 

950  READ(5, 1001 )  DELS! 1 ), ARAT 10 , EL  AS , DGAMAD 
WRITE(6 , 1100 )  DELS (l),ARATIO, EL AS, DGAMAD 
1100  FORMAT !/,3X,"DELS!l)=",E13.6, 

A3X,"ARATIO=",E13.6,3X,"ELAS=",E13.6,3X,"DGAMAD=",E13.6) 

READ( 5 ,  1001)  RB,EJ,RJ,RN 
WRITE!6,1003)RB,EJ,RJ,RN 

1003  FORMAT!/, 3X, 3HRB= , E13.6,1X,3HEJ=,E13.6,1X,3HRJ=,E13.6,1X, 3HRN= 
4E13.6) 

READC5,  1001)  RC,RO,XFA,XFF 
URITE(6,1020)RC,RO, XFA , XFF 

1020  FORMAT !/,3X,3HRC=,E13.6,lX,3HR0=,E13.6, IX, 4HXFA=, 
t  El  3 . 6  , 1X,4HXFF= , E13 . 6 ) 

READ! 5,  1001)  XLP,XCG,SN,THETND 
WRITE! 6 , 1004 )XLP ,XCG, SN , THETND 

1004  FORMAT !/,3X,4HXLP=,E13.6» IX, 4HXCG=, E13 . 6 , IX, 
*3HSN=,E13.6,1X,7HTHETND=,E13.6) 

READ! 5, 1001)  THETFD, VEL , FREKOT , DALFAD 
WRITE! 6, 1005)  THETFD, VEL , FREKOT , DAL  FAD 

1005  FORMAT!/, 3X, "THETFD=", E13.6,1X»"VEL="»E13.6» IX, "FREKOT ="> 

AEl 3 . 6 , 3X, "DALFAD=", E13. 6 ) 

READ! 5 , 1001 )  GAMSTP , YSTOP ,ZSTOP  »ZSHIFT 
READ! 5, 1001)  FRAD, FANG, DTMIN , DTMAX 
READ! 5, 1001)  WEIGHT. XICG.ALPHAD, ALSTP 
READ! 5, 1001 )  TIMEI , TIMEF, FREQI 
READ! 5, 1001)  PHIMIN, VELF.W1I .GAMMAD 
WRITE! 6 , 1006 )  GAMSTP, YSTOP, ZSTOP.ZSHIFT 
WRITE!6,1007)  FRAD , FANG, DTMIN , DTMAX 
WRITE!6,1009)  WEIGHT , XICG, ALPHAD, ALSTP 
WRITE!6,1010)  TIMEI ,TIMEF, FREQI 
WRITE!6,1011)  PHIMIN, VEL F,W1 I, GAMMAD 

1006  FORMAT!/ , 3X, 7HGAMSTP= , E13.6,2X,6HYST0P=,E13.6,2X, 
*6HZSTOP=,E13.6,2X,"ZSHIFT=",E13.6) 

1007  FORMAT!/ , 3X, "FRAD=", E13 . 6 » 2X, "FANG="» 

t El  3 . 6 , 2X, "DTMIN =" , E13 . 6 , 2X, "DTMAX ="> E13 . 6  ) 

1009  FORMAT!/ , 3X, "WEIGHT5" , El 3 . 6 , 2X, "XICG=", £13.6, 2X, "ALPHAD5" , 
*E13.6,2X,"ALSTP=",E13.6) 

1010  FORMAT !/,3X,"TIMEI=",E13.6, 2X,"TIMEF="» 

*E13.6,2X,"FREQI=",E13.6) 

1011  FORMAT !/3X,"PHIMIN=",E13.6,"  DEG.",2X,5HVELF=,E13.6, 

t  2X, 4HW1 1  =  , E13 . 6 , 2X, 7HGAMMAD= , El  3 . 6 ) 

READ!5,1002)!SNUM! I ) , DENSITY! I ) , YIELD! I ) ,ZM! I ) ,XIRD! I ) , 1=2, NLAY) 
WRITE!6 , 1301 ) 

1301  FORMAT!//, 5X, "LAYER", 12X, "SNUM", 14X, "DENSITY", 7X, 
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*"YIELD",10X,"ZM",16X,"XIRD",/,6X,"NO.",7X, 

*"< YOUNG'S  S-NUMBER)",2X,"(LB-SECt2/INt4)",3X, 

*"<PSI)",3X,"( LAYER  DEPTH-IN .  )", 2X, "(TENS . RIG. INDEX)",/, 
*/,7X,"l  *x*AIR***NO  RESISTANCE**  *") 
MR1TE(6 , 1302) (DDUM, I , SNUM( I > , DENSITY! I ),YIELD(I), 
t  ZM(I),XIRD(I),I=2»NLAY) 

1302  FORMAT(/, 4X, A1 ,I3,12X,F6.2,10X,1PE12.4,2X,E11.3,E12.3, 
t  7X,0PF9.1) 

PHIMIN=PHIMIN* DEGRAD 
SYNEFEE=SIN(PHIMIN) 

TANFEE=TAN(PHIMIN) 

XIINV-AGRAV/XICG 
W=WEIGHT;  X=0. 

XM  =  AGRAV  /  W 

DO  235  J=2,NLAY 

SHRGD(J)=XIRD(J)**<l./3.) 

IFCSNUM(J))  994,991,992 

991  LMEDIA(J)=.F. 

GO  TO  235 

992  LMEDIA(J)=.T. 

235  CONTINUE 

ZCUT=-MAX(XLP-XCG»RO*<SHRGD<NLAY)*l ) )-RO 
ZCUT  =ZCUT-ZM<  NL AY-1 ) 

IF(IOPSHP.NE.O)  GO  TO  14 
READ(5,1000)  NEL 
WRITE(6,1150)  NEL 

1150  FORMAT (// , 3X , "NEL  (NO.  OF  ELEMENTS)  =",I4) 

WRITE(6,1151) 

1151  FORMATS, 3X, "THETA-ARRAY. . .") 

READ( 5  ,1001)  (THETAD(I),I=1,NEL) 

WRITE<6,1001)  <THETAD(I),I=1,NEL) 

WRITE(6 ,1152) 

1152  FORMAT (/, 3X, "R- ARRAY .  . .") 

READC5  ,1001)  (  R( I ) , 1=1 , NEL ) 

WRITEC  6 , 1001 )  <  R( I ) ,1  =  1 , NEL ) 

URITE(6 , 1153) 

1153  FORMAT (/, 3X, "DEL S-ARR AY. .  .") 

READ( 5  ,1001)  C  DELS( I ) , 1=1 , NEL ) 

WRITE(6,1001)  (  DELSC I ), 1=1, NEL) 

14  CONTINUE 
RETURN 

994  WRITE(6,995) 

995  F0RMAT(/,2X, "NEGATIVE  S-NUMBER. . .RUN  ABORTED.") 

999  CALL  CALCMP(0.,0. ,9999,2) 

CALL  DETACH(5, , ) 

STOP 

END 


CM****  CONVERT  CALCOMP  CRT  PLOT  CALL  TO  CALCOMP  DRUM  PLOT  CALL  *** 
SUBROUTINE  CALCMP(X,Y,IZ,I) 

DATA  J0/-1/ 

IF( I . EQ. 2)  GO  TO  30 
IF(IZ. LT.2)  GO  TO  10 
CALL  PLOT(X,Y,2) 


G27 
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60  TO  20 

10  CALL  PLOT (X#Y<  3) 

60  TO  20 

30  IFCIZ.EQ.9999)  60  TO  40 
J0=-J0 

CALL  PLOT(-2.,-2.5,-3> 

CALL  PLOT ( 0 . >  11 . ,2) 

CALL  PL0T(8.5# 11 . >2) 

CALL  PLOT(8.5,0. ,2) 

CALL  PLOT (  0  .  <  0  .  #  2 ) 

IF( JO . EQ.-l )  60  TO  22 
CALL  PL0T(2. »14.5*-3) 

60  TO  20 

22  CALL  PLOT (12. >-9.5, -3) 

GO  TO  20 

40  CALL  PLOT (10. >0. >3) 

CALL  PLOT ( 0 . > 0 . , 999) 

20  RETURN 


APPENDIX  H 


NOTATION 

Differential  area  on  projectile  surface 

Net  force  components  in  the  projectile-fixed  y-  (lateral) 
and  z-  (axial)  directions,  respectively 

Net  force  components  in  the  target-fixed  Y-  and  Z- 
directions,  respectively 

Projectile's  transverse  mass  moment  of  inertia  about  its 
CG 

Free-surface  parameter 

Logic  flag  used  to  signify  that  the  end  of  the  problem  is 
approaching 

Projectile  nose  length 
Projectile  mass 

Total  moment  exerted  on  the  projectile  about  its  CG 

Unit  vector  normal  to  projectile's  surface 

Young's  nose  performance  coefficient  (Reference  2) 

Number  of  integration  time  steps  required  to  pass  through 
the  thinnest  target  layer 

Damage  radius 
Projectile  aftbody  radius 

Local  cylindrical  radius  at  a  given  point  on  the  projec¬ 
tile  surface 

Path  length  traveled  by  projectile  up  to  time  t 

Young' 8  S-number  (Reference  2),  soil  penetrability  index 

Time 

Maximum  time  allowed  for  current  problem 
Initial  time 

Thickness  of  the  thinnest  target  layer  in  a  multilayer 
problem 

Velocity  of  any  point  on  projectile  surface 
Outward  normal  component  of  local  velocity 


Projectile  CG  velocity 

Minimum  CG  velocity  allowed  for  current  problem 
Projectile  impact  velocity 
Projectile  weight 

Distance  from  projectile  nose  tip  to  CG 

Lateral  and  axial  position  in  projectile,  respectively 
Maximum  absolute  value  of  Y  allowed  for  current  problem 

CG  displacement  components  in  the  target-fixed  Y-  and  Z- 
directions 

Initial  value  of  horizontal  position  and  velocity  of  CG 

Velocity  components  of  CG  in  the  target-fixed  Y-  and  Z- 
directions 

Net  acceleration  components  in  the  target-fixed  Y-  and  Z- 
directions,  respectively 

Average  values  of  Y  and  Z  during  current  time  step 

Values  of  Y  and  Z  at  the  beginning  of  current  time 
step 

Total  path  length  traveled  by  a  penetrator  from  impact  to 
rest 

Initial  values  of  vertical  position  and  velocity  of  CG 

Vertical  shift  in  initial  projectile  position,  positive 
downwards 

Minimum  value  of  Z  allowed  for  current  problem 

Actual  integration  time-step  size 

Time  step  calculated  using  translational  parameters 

Time  step  calculated  using  rotational  parameters 

Change  in  CG  velocity  during  time  step  At 

Change  in  CG  displacement  components  in  the  target-fixed 
Y-  and  Z-directions  during  time  step  At 

Change  in  CG  velocity  components  in  the  target-fixed  Y- 
and  Z-directions  during  time  step  At 

Change  in  projectile  orientation  during  time  step  At 

Change  in  angular  velocity  during  time  step  At 


Projectile  angle  of  attack  (l.e.,  angle  between  projec¬ 
tile  axis  and  CG  velocity  vector,  positive  when  velocity 
vector  is  clockwise  from  axis) 

Initial  value  of  angle  of  attack  a 

Maximum  absolute  value  of  a  allowed  for  current  problem 

Obliquity  angle,  measured  counterclockwise  from  target- 
fixed  Z-axis  to  projectile-fixed  z-axis 

Initial  values  of  obliquity  and  angular  velocity 
Maximum  value  of  y  allowed  for  current  problem 
Angular  velocity 

Average  value  of  y  during  current  time  step 
Value  of  y  at  the  beginning  of  current  time  step 
Angular  acceleration 

Angle  between  local  tangential  velocity  component  and 
local  total  velocity  vector 

Mass  density  of  either  a  hard  target  material  or  a 
projectile  structure  material 

Local  compressive  normal  stress  on  projectile  surface 
Smoothed  normal  stress  function  (Appendix  E) 


DISTRIBUTION  LIST 


Address 


POD 

Director,  Defease  Advanced  Research  Project  Agency, 

1400  Wilson  Blvd.,  Arlington,  Virginia  22209 
ATTN:  Technical  Library 

Defense  Documentation  Center,  Cameron  Station, 

Alexandria,  Virginia  22314 
ATTN:  TC/Myer  B.  Kahn 

Director,  Defense  Intelligence  Agency,  Washington, 

D.  C.  20301 

ATTN:  Technical  Library 

Defense  Mapping  Agency,  Washington,  D.  C.  203C5 
ATTN:  H.  Lindsey 

Director,  Defense  Nuclear  Agency,  Washington,  D.  C.  20305 
ATTN:  STS I  (Archives) 

STTL  (Technical  Library) 

SPSS 

DDST 

TISI 

Under  Secretary  of  Defense  for  Research  and  Engineering, 
Department  of  Defense,  Washington,  D.  C.  20301 
ATTN:  Deputy  Director,  Tactical  Warfare 

Commander,  Field  Command,  Defense  Nuclear  Agency,  Klrtland 
Air  Force  Base,  New  Mexico  87115 
ATTN:  FCTMD 

Interservice  Nuclear  Weapons  School,  Kirtland  Air  Force 
Base,  New  Mexico  87115 
ATTN:  Technical  Library 

Director,  Joint  Strategic  Target  Planning  Staff,  JCS, 
Offutt  Air  Force  Base,  Nebraska  68113 
ATTN:  Science  and  Technology  Info  Library 


Assistant  Chief  of  Staff  for  Force  Development,  Department 
of  the  Army,  Washington,  D.  C.  20310 
ATTN:  Technical  Library 

Commander,  U.  S.  Army  Corps  of  Engineers,  Washington, 

D.  C.  20314 
ATTN:  DAEN-ASI-L 
DAEN-MPE-T 
DAEN-RDL 
DAEN-RDM 

Chief  of  Staff  for  Research,  Development,  and  Acquisition, 
Department  of  the  Army,  Washington,  D.  C.  20310 
ATTN:  Technical  Library 

Commander,  Harry  Diamond  Laboratories,  Department  of  the 
Army,  2800  Powder  Mill  Road,  Adelphi, .Maryland  20783 
ATTN:  AMXDO-RBH/J.  Gwaltney 

Commanding  Officer,  Picatinny  Arsenal,  Dover, 

New  Jersey  07801 
ATTN:  Technical  Library 

Director,  U.  S.  Army  Ballistic  Research  Laboratories, 
Aberdeen  Proving  Ground,  Maryland  21005 
ATTN:  G.  Roecker 

Technical  Library /Edward  Baicy 

Division  Engineer,  U.  S.  Army  Engineer  Division, 

Huntsville,  P.  0.  Box  1600,  West  Station,  Huntsville, 
Alabama  35807 

ATTN:  HNDED-R/Michael  M.  Dembo 

Commander,  Army  Materials  and  Mechanics  Research  Center, 
Watertown,  Massachusetts  02172 
ATTN :  Technical  Library 

Commander,  U.  S.  Army  Material  Development  and  Readiness 
Command,  5001  Elsenhower  Avenue,  Alexandria, 

Virginia  22333 
ATTN :  Technical  Library 

Commander,  U.  S.  Army  Missile  Command,  Redstone  Arsenal, 
Alabama  35809 
ATTN:  Technical  Library 

DRCPM-PE-X  (William  K.  Jann) 


Army  (Continued) 

Commander,  U.  S.  Army  Nuclear  and  Chemical  Agency, 

7500  Backlick  Road,  Bldg  2073,  Springfield, 

Virginia  22150 

ATTN:  Technical  Library  1 

Commander,  U.  S.  Army  Armament  Command,  Rock  Island 
Arsenal,  Rock  Island,  Illinois  61201 
ATTN:  Technical  Library  1 


Navy 


Chief  of  Naval  Research,  Department  of  the  Navy, 

Arlington,  Virginia  22217 

ATTN:  Technical  Library  1 

Naval  Civil  Engineering  Laboratory,  Port  Hueneme, 

California  93043 

ATTN:  Technical  Library  1 

Naval  Facilities  Engineering  Command,  200  Stovall  Street, 

Alexandria,  Virginia  22332 

ATTN:  Technical  Library  1 

Director,  Naval  Research  Laboratory,  Washington, 

D.  C.  20375 

ATTN:  Code  2027 /Technical  Library  1 

Commander,  Naval  Surface  Weapons  Center,  White  Oak,  Silver 
Spring,  Maryland  20910 

ATTN:  Technical  Library,  Code  730  1 

Commander,  Naval  Surface  Weapons  Center,  Dahlgren 
Laboratory,  Dahlgren,  Virginia  22448 

ATTN:  Technical  Library  1 

Commander,  Naval  Weapons  Center,  China  Lake, 

California  93555 

ATTN:  Code  2663,  Dr.  Carl  Austin  1 

Commander,  Naval  Weapons  Evaluation  Facility,  Kirtland  Air 
Force  Base,  Albuquerque,  New  Mexico  87117 

ATTN:  Technical  Library  1 


ill 


Address 


No.  of 
Copies 


Air  Force 


Air  Force  Armament  Laboratory  (AFSC) ,  Eglin  Air  Force 
Base,  Florida  32542 
ATTN:  DLYV  (Mr.  John  Collins) 

DLJW  (Mr.  G.  R.  Griner) 

Technical  Library 

Air  Force  Weapons  Laboratory  (AFSC) ,  Klrtland  Air  Force 
Base,  New  Mexico  87117 
ATTN :  Technical  Library 

Headquarters,  Air  Force  Systems  Command,  Andrews  Air  Force 
Base,  Washington,  D.  C.  20331 
ATTN :  Technical  Library 

Oklahoma  State  University  Field  Office,  P.  0.  Box  1925, 
Eglin  Air  Force  Base,  Florida  32542 
ATTN:  Ed  Jackett 

Air  Force  Institute  of  Technology,  Air  University, 
Wright-Patterson  Air  Force  Base,  Ohio  45433 
ATTN :  Library 


Other  Government  Agencies 

Lawrence  Livermore  National  Laboratory,  P.  0.  Box  808, 
Livermore,  California  94550 
ATTN:  Technical  Library 

Nuclear  Regulatory  Commission,  Directorate  of  Licensing 
Regulations,  Washington,  D.  C.  20555 
ATTN:  Mr.  Lawrence  Shao,  Chief,  Struc  Engrg  Br 

Sandia  National  Laboratories,  P.  0.  Box  5800,  Albuquerque, 
New  Mexico  87115 
ATTN:  Mr.  William  Patterson 
C.  Wayne  Young 


POD  Contractors 


Avco  Corporation,  201  Lowell  Street,  Wilmington, 
Massachusetts  01887 

ATTN:  David  Henderson,  Government  Products  Group 
Research  Library  A830,  Room  7201 


•  •  •  9 


i 

i 

i 


i 


i 


i 


i 


i 


l 


i 

i 


i 

i 


iv 


Address 


No.  of 
Copies 


1 

c* 


1 


POD  Contractors  (Continued) 

Boeing  Aerospace  Company,  F.  0.  Box  3999,  Mall  Stop  13-11, 
Seattle,  Washington  98124 
ATTN:  Aerospace  Library 

Braddock,  Dunn,  &  McDonald,  Inc.,  1920  Aline  Avenue, 
Vienna,  Virginia  22180 
ATTN :  Technical  Library 

California  Research  and  Technology,  Inc.,  6269  Variel 
Avenue,  Woodland  Hills,  California  91367 
ATTN:  Dr.  K.  N.  Kreyenhagen 

Institute  for  Defense  Analyses,  400  Army-Navy  Drive, 
Arlington,  Virginia  22202 
ATTN:  IDA  Librarian  Ruth  S.  Smith 

Lockheed  Missiles  and  Space  Company,  Inc.,  3251  Hanover 
Street,  Palo  Alto,  California  94304 
ATTN:  Tech  Info  Center  D/C0LL 

Martin  Marietta  Aerospace,  Orlando  Division,  P.  0.  Box 
5837,  Orlando,  Florida  32805 
ATTN:  M.  Anthony 

New  Mexico  Engineering  Research  Institute,  University  of 
New  Mexico,  Box  25,  University  Station,  Albuquerque, 

New  Mexico  87131 
ATTN:  G.  E.  Triandafalidis 

Pacifica  Technology,  P.  0.  Box  148,  Del  Mar, 

California  92014 
ATTN:  R.  L.  Bjork 

Physics  International  Company,  2700  Merced  Street,  San 
Leandro,  California  94577 
ATTN:  Doc  Control  for  Mr.  Larry  Behrman 

R&D  Associates,  P.  0.  Box  9695,  Marina  del  Rey, 

California  90291 
ATTN :  Technical  Library 

TRW  Defense  &  Space  Systems  Group,  Mail  Station  Rl/2154, 
One  Space  Park,  Redondo  Beach,  California  90278 
ATTN:  Technical  Info  Center /S-1930 


1 


1 


1 


1 


1 


1 


1 


1 


1 


1 


1 


v 


In  accordance  with  letter  fro*  DAEN-RDC,  DAEN-ASI  dated 
22  July  1977,  Subject:  Facsiaile  Catalog  Cards  for 
Laboratory  Technical  Publications,  a  facsiaile  catalog 
card  in  Library  of  Congress  MARC  format  is  reproduced 
below. 
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"Prepared  for  Defense  Nuclear  Agency  under  DNA 
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